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HIGH-ORDER/SPECTRAL  METHODS  ON  UNSTRUCTURED  GRIDS 
I.  TIME-DOMAIN  SOLUTION  OF  MAXWELL’S  EQUATIONS  * 

J.S.  HESTHAVEN'i  AND  T.  WARBURTON+ 


Abstract.  We  present  an  ab  initio  development  of  a  convergent  high-order  accurate  scheme  for  the 
solution  of  linear  conservation  laws  in  geometrically  complex  domains.  As  our  main  example  we  present  a 
detailed  development  and  analysis  of  a  scheme  suitable  for  the  time-domain  solution  of  Maxwell’s  equations 
in  a  three-dimensional  domain.  The  fully  unstructured  spatial  discretization  is  made  possible  by  the  use  of  a 
high-order  nodal  basis,  employing  multivariate  Lagrange  polynomials  defined  on  the  triangles  and  tetrahedra. 
Careful  choices  of  the  unstructured  nodal  grid  points  ensure  high-order/spectral  accuracy,  while  the  equations 
themselves  are  satisfied  in  a  discontinuous  Galerkin  form  with  the  boundary  conditions  being  enforced  weakly 
through  a  penalty  term.  Accuracy,  stability,  and  convergence  of  the  semi-discrete  approximation  to  Maxwell’s 
equations  is  established  rigorously  and  bounds  on  the  global  divergence  error  are  provided.  Concerns  related 
to  efficient  implementations  are  discussed  in  detail. 

This  sets  the  stage  for  the  presentation  of  examples,  verifying  the  theoretical  results,  as  well  as  illustrating 
the  versatility,  flexibility,  and  robustness  when  solving  two-  and  three-dimensional  benchmarks  in  computa¬ 
tional  electromagnetics.  Pure  scattering  as  well  as  penetration  is  discussed  and  high  parallel  performance  of 
the  scheme  is  demonstrated. 

Subject  classification.  Applied  Mathematics 

Key  words,  high-order /spectral  accuracy,  stability,  convergence,  unstructured  grids,  Maxwell’s  equa¬ 
tions 

1.  Introduction.  The  ability  to  accurately  and  reliably  model  wave-dominated  problems  continues 
to  be  an  essential,  and  in  many  cases  an  enabling,  technology  in  the  development  and  analysis  of  emerging 
technologies  such  as  stealth  technology,  noise  reduction,  subsurface  exploration  and  optical  communication  to 
name  a  few.  These  are  all  problems  characterized  by  being  very  large  in  terms  of  a  characteristic  wavelength, 
geometrically  extremely  complex,  often  composed  of  a  heterogeneous  collection  of  different  materials  and  all 
requiring  a  high  fidelity  solution  with  a  rigorous  control  of  the  numerical  errors.  Even  for  linear  problems 
such  conditions  forces  one  to  look  beyond  standard  computational  techniques  and  seek  new  computational 
frameworks  enabling  the  accurate,  efficient,  and  robust  modeling  of  wave-phenomena  over  long  times  in 
settings  of  a  realistic  geometric  complexity. 

The  requirement  that  one  can  accurately  propagate  waves  over  many  periods  of  time  naturally  suggests 
that  high-order /spectral  methods  be  considered  [1].  On  the  other  hand,  the  use  of  such  methods  is  tra¬ 
ditionally  in  conflict  with  the  need  for  significant  geometric  flexibility  by  being  restricted  to  fairly  simple 
geometries.  The  standard  approach  to  overcome  this  restriction  is  to  introduce  a  multi-element  formulation 
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in  which  the  basic  building  block  is  parametrically  mapped  cubes  in  the  spirit  of  finite  element  methods. 
This  approach  has  been  very  successfully  applied  to  the  solution  of  problems  in  fluid  mechanics  [2,  3,  4], 
gasdynamics  [5,  6,  7,  8,  9,  10],  and  electromagnetics  [11,  12,  13,  14,  15]. 

While  such  techniques,  when  applicable,  are  powerful  they  do  suffer  from  the  need  to  tile  the  computa¬ 
tional  using  only  hexahedral  elements.  Unfortunately,  automated  grid  generation  using  only  such  elements 
for  general  three-dimensional  computational  problems  of  a  realistic  complexity  remains  a  very  nontrivial 
task  and  is  typically  very  time-consuming.  Furthermore,  spatial  adaptation,  while  certainly  possible,  is 
quite  a  challenge  with  a  method  based  solely  on  hexahedral  elements.  On  the  other  hand,  automated  grid 
generation  employing  a  fully  unstructured  grid  is  significantly  more  mature,  due  mainly  to  extensive  devel¬ 
opments  within  the  finite-element  community.  Spatial  grid  adaptation  is  also  considerably  easier  within  a 
fully  unstructured  grid  formulation. 

It  is  with  these  issues  in  mind  that  we  present  an  ab  initio  development  of  a  computational  framework 
that  combines  the  strengths  of  a  high-order/spectral  formulation  with  the  flexibility  of  a  fully  unstructured 
grid.  The  formulation  relies  on  the  resolution  of  two  central  issues.  On  one  hand  we  shall  discuss  in  detail 
how  to  represent  functions  defined  on  triangles  and  tetrahedra  to  high  accuracy  and  how  this  translates  into 
the  construction  of  basic  operators  needed  to  solve  partial  differential  equations.  On  the  other  hand  we  need 
to  address  the  issue  of  how  to  use  such  a  high-order  representation  to  formulate  a  convergent  scheme  suitable 
for  solving  systems  of  linear  hyperbolic  problems  in  general  and  Maxwell’s  equations  in  particular. 

Much  in  the  spirit  of  the  original  work  on  spectral  element  methods  [2,  3]  we  shall  focus  on  the  formulation 
of  efficient  and  flexible  unstructured  grid  methods  using  nodal  elements.  This  is  in  contrast  to  past  attempts 
to  develop  high-order  unstructured  grid  methods,  suitable  for  solving  time  dependent  problems,  which  have 
been  focused  on  the  use  of  high-order  modal  expansions,  e.g.,  [16,  17,  18,  19,  20,  21].  In  these  works,  modal 
expansions  of  orthogonal  polynomials  defined  on  the  simplex  are  utilized  while  a  straightforward  monomial 
basis  is  used  in  [22]  (see  also  [23]  and  references  therein)  much  in  the  tradition  of  classical  high-order  finite 
element  methods  for  elliptic  problems  [24,  25]. 

In  contrast  to  the  classical  spectral  element  approach,  however,  we  do  not  seek  a  globally  continuous 
solution  but  rather  require  that  the  equations  be  satisfied  in  a  discontinuous  Galer kin/penalty  fashion.  This 
is  related  to  the  classic  discontinuous  Galerkin  finite  element  method  [23]  although  the  present  approach 
represents  a  more  general  formulation,  containing  the  classic  discontinuous  Galerkin  formulation  as  a  special 
case.  Such  more  general  techniques  have  been  known  in  the  context  of  spectral  methods  as  penalty  meth¬ 
ods  [26]  for  a  while  and  recently  stable  formulations  on  general  one-dimensional  [27],  triangular  [28],  and 
tetrahedral  domains  [29]  have  been  discussed.  These  methods  all  share  the  great  advantage  of  a  complete 
decoupling  of  all  elements,  hence  enabling  high  parallel  efficiency,  and  allows  for  discontinuous  solutions 
between  elements  in  a  natural  way.  As  we  shall  see  later,  this  is  essential  in  allowing  for  the  inclusion  of 
material  interfaces  in  a  natural  and  straightforward  manner. 

While  the  majority  of  what  we  shall  discuss  is  of  a  very  general  nature  we  have  chosen  to  discuss  in  detail 
the  development  and  analysis  of  a  high-order/spectral  accuracy  unstructured  grid  scheme  for  the  solution 
of  Maxwell’s  equations  in  the  time-domain.  This  is  not  only  a  challenging  problem  but  also  a  problem  of 
significant  contemporary  interest  due  to  emerging  technologies  such  as  broad-band  target  illumination  and 
penetration,  advanced  materials  and  diffraction  based  modern  optics,  all  characterized  by  being  electrically 
large,  having  a  significant  separation  of  scales  and  requiring  substantial  geometric  flexibility  of  the  compu¬ 
tational  framework.  On  the  other  hand,  Maxwell’s  equations  serve  as  an  excellent  example  of  numerous 
other  linear  hyperbolic  systems  of  equations  in,  e.g.,  elasticity,  acoustics,  solid  mechanics  etc,  for  which 
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the  presented  framework  can  be  adapted  with  little  effort.  In  part  II  of  this  work  [30]  we  shall  discuss  in 
detail  generalizations  of  the  proposed  computational  framework  with  an  emphasis  on  the  solution  of  general 
systems  of  conservation  laws. 

What  remains  of  the  paper  is  organized  as  follows.  In  Sec.  2  we  set  the  stage  by  briefly  describ¬ 
ing  the  physical  setting,  Maxwell’s  equations,  their  normalized  and  scattered  field  formulations,  as  well  as 
boundary  conditions  at  material  interfaces  and  metallic  boundaries.  The  first  step  in  the  construction  of 
a  high-order /spectral  unstructured  grid  scheme  for  the  solution  of  Maxwell’s  equations  is  taken  in  Sec.  3 
where  we  introduce  a  Lagrangian  high-order  basis  on  the  general  curvilinear  simplex.  In  the  appendix  we 
include  a  discussion  of  techniques  allowing  for  efficient  and  accurate  implementations  of  the  basic  operators, 
e.g.,  differentiation,  filtering,  and  high-order  integration  in  volumes  and  on  faces.  By  providing  the  basic 
building  block  for  the  spatial  approximation,  this  development  sets  the  stage  for  the  formulation  of  a  high- 
order/spectral  convergent  scheme  for  solving  Maxwell’s  equations  as  discussed  in  Sec.  4.  The  convergence 
of  the  scheme,  being  a  generalized  discontinuous  Galer kin/penalty  method,  is  established  in  the  classic  way 
through  consistency  as  well  as  local  and  global  stability.  A  stronger  and  optimal  result  is  furthermore  estab¬ 
lished  by  showing  the  scheme  to  be  error-bounded,  guaranteeing  at  most  linear  growth  in  time  and  control 
over  the  growth  rate.  This  result  is  also  used  to  establish  bounds  on  the  behavior  of  the  divergence  error. 
Verification  and  performance  of  the  complete  scheme  is  the  topic  of  Sec.  5  where  we  present  a  number  of 
simple  tests,  verifying  the  theoretical  results,  prior  to  illustrating  the  efficiency,  versatility,  and  robustness 
of  the  computational  framework  for  the  solution  of  two-  and  three-dimensional  scattering  and  penetration 
problems.  We  shall  also  briefly  discuss  measures  taken  in  the  implementation  of  the  scheme  to  ensure  efficient 
execution  on  large  scale  contemporary  parallel  computational  platforms.  In  Sec.  6  we  conclude  by  offering 
a  few  remarks  and  guidelines  for  future  work  within  the  present  framework. 

2.  The  Physical  Setting  and  Maxwell’s  Equations.  We  shall  concern  ourselves  with  the  direct 
solution  of  Maxwell’s  equations  on  differential  form 


(1) 


dD 

~dt 


=  V  xH  +  J 


dB 

~dt 


— V  x  E  , 


(2)  V  •  D  =  p  ,  V  ■  B  =  0  , 

within  the  general  three-dimensional  domain,  f l,  with  the  charge  distribution,  p(x,t).  The  electric  field, 
E(x,t ),  and  the  electric  flux  density,  D(x,t),  as  well  as  the  magnetic  field,  H(x,t),  and  the  magnetic  flux 
density,  B(x,t),  are  related  through  the  constitutive  relations 

D  =  sE  ,  B  =  (1H  . 

The  permittivity  tensor,  e,  as  well  as  the  permeability  tensor,  fi.  are  in  general  anisotropic  and  may  depend 
on  space  and  time  as  well  as  the  strength  of  the  fields  themselves.  The  current,  J,  is  typically  assumed  to 
be  related  to  the  electric  field,  E.  through  Ohms  law,  J  =  aE,  where  d  measures  the  finite  conductivity, 
although  more  complex  relations  are  possible. 

In  this  work,  we  shall  restrict  the  attention  to  materials  which  can  be  assumed  isotropic,  linear  and 
time-invariant,  in  which  case  the  constitutive  relations  take  the  form 


D  —  SQSrE  ,  B  —  fi.0prH  . 


3 


Here  io  =  8.854  x  10-12  F/m  and  po  =47 rx  10-7  H/m  represent  the  vacuum  permittivity  and  permeability, 
respectively,  and  er(x)  and  y ir(x )  refers  to  the  relative  permittivity  and  permeability,  respectively,  of  the 
materials. 

Taking  the  divergence  of  Eq.(l)  and  applying  Eq.(2)  in  combination  with  Gauss’  law  for  charge  conser¬ 
vation  immediately  confirms  that  if  the  initial  conditions  satisfy  Eq.(2),  and  the  fields  are  evolved  according 
to  Maxwell’s  equations,  Eq.(l),  the  solution  will  satisfy  Eq.(2)  at  all  times.  Hence,  one  can  view  Eq.(2)  as  a 
consistency  condition  on  the  initial  conditions  and  limit  the  solution  to  the  time-dependent  part  of  Maxwell’s 
equations,  Eq.(l). 

To  simplify  matters  further,  we  shall  consider  the  non-dimensionalized  equations  for  which  we  introduce 
the  normalized  quantities 


t 

t  =  -r - 

L/c0 


where  L  is  a  reference  length,  and  Co  =  (§oPo)  1 represents  the  dimensional  vacuum  speed  of  light.  The 
fields  themselves  are  normalized  as 

E=ZX°  ,  H  =  t  ,  J=  , 

H0  H0  H0 


where  Z0  =  \J po /so  refers  to  the  dimensional  free  space  intrinsic  impedance,  and  Ho  is  a  dimensional 
reference  magnetic  field  strength. 

With  this  normalization  Eq.(l)  takes  the  nondimensional  form 


(3)  £r~Qt  =  V  x  H  +  J  ’  tJtr~dT  =  _V  X  E  ’ 

which  is  the  general  form  of  the  equations  we  consider  in  the  following. 

To  solve  Maxwell’s  equations  in  the  vicinity  of  boundaries,  penetrable  or  not,  we  shall  need  boundary 
conditions  relating  the  field  components  on  either  side  of  the  boundary. 

Assuming  that  a  normal  unit  vector,  n,  to  the  boundary  is  given,  the  boundary  conditions  on  the  electric 
field  components  take  the  form 


n  x  [Ei  -  E2)  =  0:  ,  n  {D x-  D2)  =  ps  , 

where  Et  and  D, ,  i  =  (1,2),  represent  the  fields  on  either  side  of  the  interface  and  ps  represents  a  surface 
charge.  Equivalently,  the  conditions  on  the  magnetic  fields  are  given  as 

n  x  (Hi  -  H2)  =  Js  ,  n-(B1-  B2)  =  0  , 

where  Js  represents  a  surface  current  density. 

In  the  general  case  of  materials  with  finite  conductivity,  no  surface  charges  and  currents  can  exist,  and 
the  simplified  conditions  take  the  form 


(4)  n  x  (Ei  -  E2)  =  0  ,  nx(Hi-H 2)  =  0  , 

expressing  continuity  of  the  tangential  field  components,  while  the  normal  components  of  the  flux  densities 
must  satisfy 
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(5) 


ii  ■  (D i  -  D 2)  =  0  ,  ii  ■  (B1  -  B2 )  =  0  , 


i.e.,  they  are  likewise  continuous,  while  the  normal  components  of  the  fields  themselves  are  discontinuous. 

For  the  important  special  case  of  a  perfect  conductor,  the  conditions  take  a  special  form  as  the  perfect 
conductor  supports  surface  charges  and  currents  while  the  fields  are  unable  to  penetrate  into  the  body,  i.e., 


(6)  nxE  =  0,n-B=0. 

2.1.  The  Scattered  Field  Formulation.  For  scattering  and  penetration  problems  involving  linear 
materials  it  is  often  advantageous  to  exploit  the  linearity  of  Maxwell’s  equations  and  solve  for  the  scattered 
field,  ( ES,HS ),  rather  than  for  the  total  field,  ( E,H ),  which  are  trivially  related  as 

E  =  E7  +  Es  ,  H  =  H'  +  Hs  , 

where  {E\H2)  represents  the  incident  field,  illuminating  the  scattering  object.  Assuming  that  (E\H2) 
represents  a  particular  solution  to  Maxwell’s  equations,  one  recovers  the  scattered  field  formulation 

(7)  e, =  VxHs+aEs-  (er  -  e* )  +  (a  -  a^E*  , 

f)  FTS  f)  FTi 

(8)  =-vxF-(fr-/i’)—  , 

where  e’r{x),  n7r(x),  and  cr7(x)  refers  to  the  relative  permittivity,  permeability  and  conductivity  of  the  media 
in  which  the  incident  field  represents  a  solution  to  Maxwell’s  equations.  To  simplify  matters  we  have  assumed 
Ohms  law,  J  =  aE.  We  note  that  the  important  special  case  of  a  vacuum  field  illuminating  the  scattering 
object  is  recovered  by  using  e*.  =  /j’r  =  1,  a7  =  0,  and  using  a  free  space  solution  in  the  forcing  function. 

In  this  formulation,  the  boundary  conditions  along  a  dielectric  interface  take  the  form 


(9)  n  x  {E{  —  E2)  =  0  ,  n  x  (H\  —  Hs2)  =  0  , 

for  the  tangential  components,  while  the  conditions  on  the  scattered  field  components  becomes 


(10)  n  x  Es  =  —n  x  E7  ,  n  ■  Bs  =  -pTn  •  H7  , 

in  the  case  of  a  perfectly  conducting  boundary.  As  we  shall  see  shortly,  there  is  no  need  to  consider  the 
conditions  on  the  normal  components  further. 

3.  The  Nodal  Element.  We  shall  seek  approximate  solutions  to  Maxwell’s  equations  in  a  general 
domain,  f l,  possibly  containing  a  heterogeneous  collection  of  scattering  and  penetrable  bodies.  To  facilitate 
the  required  geometric  flexibility,  we  represent  the  computational  domain  as  the  union  of  K  non-overlapping 
body-conforming  d-simplices,  D.  Hence,  for  two-dimensional  problems  we  shall  use  triangles  as  the  geometric 
building  block  while  the  tetrahedron  is  employed  to  fill  the  computational  volume. 
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Fig.  1.  Mapping  between  the  curvilinear  tetrahedral  D,  and  the  standard  tetrahedral  I,  including  the  numbering  and 
notation  employed  in  the  text. 

While  this  multi  element  formulation  is  essential  in  enabling  the  solution  of  geometrically  complex 
problems,  it  also  introduces  new  complications,  the  understanding  and  resolution  of  which  are  at  the  heart 
of  the  construction  of  the  scheme.  In  particular,  the  use  of  simplices  requires  an  understanding  of  how  to 
construct  high-order  accurate  Lagrange  interpolation  polynomials  on  such  elements  and,  subsequently,  how 
we  can  formulate  approximations  to  basic  operations  such  as  interpolation,  differentiation  and  integration 
of  functions  defined  on  general  curvilinear  d-simplices.  These  are  issues  we  shall  deal  with  in  the  following. 
For  continuity  we  shall  postpone  the  discussion  of  practical,  yet  essential,  techniques  for  the  efficient  and 
accurate  implementation  of  the  basic  operations  to  the  appendix. 

The  equally  important  question  of  how  to  exploit  this  knowledge  to  construct  global  high-order / spectral 
accuracy  solution  techniques  suitable  for  Maxwell’s  equations  as  well  as  other  linear  hyperbolic  systems  is 
the  central  issue  addressed  in  Section  4. 

3.1.  The  Curvilinear  d-Simplex.  We  start  by  assuming  that  the  computational  domain,  0,  is  de¬ 
composed  into  curvilinear  d-simplices,  D  C  Rd,  as  illustrated  in  Fig.  1  by  a  3-simplex,  a  tetrahedron.  For 
generality  we  shall  limit  much  of  the  discussions  to  the  three-dimensional  case  and  regard  the  two-dimensional 
problem  as  a  natural  simplification. 

While  we  shall  not  require  that  the  faces  of  the  tetrahedron  are  planar,  such  an  assumption  will,  as  we 
shall  see  shortly,  significantly  simplify  matters  in  terms  of  analysis  as  well  as  implementation.  It  should  also 
be  noted  that  for  most  computational  problems,  the  vast  majority  of  the  elements  will  have  planar  faces 
which  thus  supplies  the  single  most  important  special  case. 

Let  us  introduce  the  standard  tetrahedron,  I  C  R3,  given  by  the  vertices 


'  -1  ' 

1 

'  -1  ' 

'  -1  ' 

Vi  = 

-1 

,  Vu  = 

-1 

,  Will  = 

1 

,  Wiv  = 

-1 

_  -1  _ 

_  -1  _ 

_  -1  _ 

1 

as  illustrated  in  Fig.  1  with  the  corresponding  vertices  in  D  termed  V1-V4.  To  fix  the  notation  within  the 
tetrahedron,  let  us  also  name  the  face  in  D  opposite  vertex  Vi,  i.e,  spanned  by  the  three  vertices  Vo,  W3,  and 
Vi,  for  face  ’a’,  that  opposite  of  vertex  v-2  for  face  ’b’  and  so  forth.  In  general  we  shall  name  the  coordinates 
in  the  physical  simplex,  D,  as  x  =  ( x,y,z )  while  the  coordinates,  £  C  I.  shall  be  referred  to  as  £  = 
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To  relate  operations  on  D  to  those  on  I  we  need  to  construct  a  smooth  and  invertible  mapping,  :  D  >  I . 
that  uniquely  relates  the  two  simplices  as  illustrated  in  Fig.  1.  In  the  case  of  a  general  curvilinear  mapping, 
this  can  be  constructed  directly  by  the  use  of  linear  transfinite  blending  functions.  Although  lengthy, 
expressions  of  these  mappings  are  straightforwardly  arrived  at  by  blending  parameterized  versions  of  faces, 
edges,  and  the  vertex-coordinates.  For  a  detailed  account  of  this  we  refer  to  [21]. 

A  particularly  important  and  simple  case  is  that  of  D  being  straightfaced  in  which  case  the  mapping 
becomes 


(11) 


*  =  ^(C)  =  - 


1  +  C  +  r]  +  C 


Vi  + 


i  +  C 


V-2  + 


1  +  1] 


V3  + 


1  +  C 


-V,i 


derived  directly  by  exploiting  that  any  point  in  the  straightfaced  tetrahedron  can  be  expressed  as  a  convex 
sum  of  the  vertices  with  the  weights  being  the  barycentric  coordinates  (see  e.g.  [21]). 

Once  the  mapping,  has  been  established,  we  can  utilize  this  to  compute  the  curvilinear  metric  of 

the  transformation  by 


dx  d£  _  <94>(C)  d£ 
<9C  dx  <9£  dx 


;rc 

xn 

xc 

vt 

Vr, 

VC 

zn 

zc  . 

C.r  Cy  C~ 

>lx  Vy  Cz 


'  1 

0 

o 

= 

0 

1 

0 

o 

0 

1  _ 

Within  this  new  metric,  the  divergence  of  a  vector  field,  F  =  (Fx,Fy,Fz),  is  expressed  on  the  well  known 
form 


V -F  =  4 


J  [dC  dr 7 

where  we  have  introduced  the  transformation  Jacobian 

I  dx 


l(JF  ■  vc)  +  !:( JF  ■  VO  +  ^(  JF  •  VC) 


J  = 


d£ 


vc  •  (V/?  x  vc) 


The  metric  also  immediately  gives  outward  pointing  normal  vectors  at  the  4  faces  of  D  on  the  form 


na  =  VC  +  Vr)  +  VC  , 


nb  =  -VC  ,  nc  =  -Vi i  ,  nd  =  -V(  ■ 


It  is  worth  while  paying  attention  to  the  special  case  of  the  mapping  between  straightsided  tetrahedra, 
Eq.(ll),  in  which  case  we  realize  that 


dx  _  94>(C)  _  1 

5C  =  <9C  =  2 


T  T 

~V(  +1>2 

rr  t 
-v{  +v:\ 

rT  rr 

-v{  +  v\ 


is  constant.  Thus,  the  full  metric,  VC,  Vr],  and  VC,  is  constant  as  is  the  transformation  Jacobian,  J,  i.e., 
every  two  straightsided  tetrahedra  are  connected  through  a  simple  linear  transformation.  As  we  discuss  in 
detail  in  the  appendix,  this  observation  can  be  exploited  to  significantly  simplify  the  implementation  of  the 
general  unstructured  scheme  by  introducing  template  operators. 

Let  us  finally  define  a  number  of  different  inner  products  on  the  curvilinear  simplex,  D.  Consider  the 
two  smooth  functions,  /[ D]  G  C[D]  and  </[D]  €  C[D]  for  which  f(x)  :  D  — >  R  and  g(x)  :  D  — »•  R.  The  global 
inner  product,  the  associated  T2-norm  and  the  inner  product  over  the  surface  of  D  are  defined  as 
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(f,9)  D=  [  f(x)g(x)dx  ,  (/,/)D  =  II/IId  ,  (f,g)SD=(f  f(x)g(x)dx  . 

J  D  JS  D 

These  local  inner  products  and  norms  form  the  basis  for  the  corresponding  global  broken  measures  as 

(/.sOfr=  d*  >  (/,/)a  =  E II/IId^  ll/llo  * 

k  k 

U,g)m  =  E  /  f(x)g(x)dx  , 

kJSDh 

where  A'  represents  the  total  number  of  elements  used  to  cover  il. 

3.2.  A  Multivariate  Polynomial  Basis  on  the  d- Simplex.  With  the  curvilinear  framework  in  place 
we  can  now  focus  the  attention  on  the  development  of  a  high-order/spectral  representation  of  a  function 
defined  on  the  elemental  element,  I,  rather  than  a  general  D. 

Contrary  to  the  approach  taken  in  [17,  21],  where  a  purely  modal  approximation  is  utilized,  we  shall  em¬ 
ploy  a  purely  nodal  scheme.  Hence,  we  assume  that  the  unknown  solutions,  q(£,t ),  can  be  well  approximated 
as 

N 

3=0 

where  Lj(£)  is  the  genuine  three-dimensional  multivariate  Lagrange  interpolation  polynomial,  Lj(£)  €  Pf,, 
where 


Pn  =  span  {f  i,j,  k  >  0]i  +  j  +  k  <  n }  , 

based  on  the  N?r  =  N  +  1  nodal  points,  £  •,  given  in  the  interior  as  well  as  on  the  boundary  of  I.  It  is 
straightforward  to  see  that  the  minimum  number  of  nodal  points  that  will  allow  this  basis  to  be  complete  is 

A^  =  I(n  +  l)(n  +  2)(n  +  3)  , 

where  n  signifies  the  maximum  order  of  the  polynomial. 

The  crucial  choice  of  a  nodal  set,  well  suited  for  Lagrange  interpolation  within  the  tetrahedron,  is  an 
issue  that  has  received  some  attention  lately  with  such  nodal  sets  being  given  in  [31]  and  [29].  The  former 
is  derived  by  using  an  minimization  procedure  for  the  identification  of  the  nodal  set  that  minimizes  an 
approximation  to  the  Lebesque  constant  while  the  approach  taken  in  the  latter  work  involves  the  solution  of 
an  electrostatic  problem  within  the  tetrahedron.  Either  procedure  results  in  fully  unstructured  nodal  sets,  an 
example  of  which  is  given  in  Fig.  2,  with  a  large  degree  of  symmetry,  exactly  nodes  within  the  tetrahedron 
and  a  very  well  behaved  Lagrange  polynomial  as  measured  through  the  growth  of  the  associated  Lebesque 
constant.  Furthermore,  both  nodal  sets  include  the  4  vertices  in  I  and  have  exactly  i(n  +  l)(n  +  2)  nodes 
at  each  of  the  four  faces.  This  latter  property  is  important  as  it  ensures  that  a  complete  two-dimensional 
polynomial  is  supported  by  the  nodes  on  each  face. 

In  this  work  we  have  chosen  to  use  the  nodal  set  from  [29]  as  the  nodes  on  which  the  Lagrange  interpo¬ 
lation  polynomials  are  based.  These  nodal  sets  are  given  for  n  up  to  10,  corresponding  to  Nf0  =  286  nodal 
points  within  each  tetrahedron  and  66  nodal  points  at  each  face. 
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Fig.  2.  Example  of  nodal  set  for  a  5th  order  interpolation,  i.e.,  N ?  —  56  nodes  within  the  tetrahedron.  In  a)  we  show  a 
3D  view  of  the  nodes  within  the  tetrahedron  while  b )  gives  a  top  view  emphasizing  the  high  degree  of  symmetry  associated  with 
the  nodal  set. 

Once  we  have  identified  a  proper  nodal  set,  we  can  proceed  with  the  formulation  of  the  interpolation 
which  must  have  the  property 


A  f(ti)  =  f(tj)  , 

for  any  /  G  C[l] .  For  the  actual  construction  of  the  interpolation  polynomials,  let  us  introduce  the  complete 
polynomial  basis,  pi(£)  £  Pf,  and  express  the  interpolation  property  as 


(12) 


N 


Vi 


/(O  =  fjPjtei) 


3= 0 


or  in  compact  form 


V/  =  f  , 

where  f  =  [/o,..,/jv]T  is  the  vector  of  expansion  coefficients,  /  =  [/(£0),  J(^n)]T  is  the  grid  vector 

and  Vij  =  pj(£j)  is  the  multi-dimensional  generalization  of  the  Vandermonde  matrix.  Clearly,  for  the 
interpolation  to  exist,  V  must  be  nonsingular  which  is  a  property  that  depends  solely  on  the  nodal  sets.  For 
polynomial  interpolation  along  the  line  it  is  well  known  that  |V|  ^  0  provided  that  the  nodes  are  distinct. 
Unfortunately,  no  such  simple  results  are  known  for  polynomial  interpolation  in  I  and  we  shall  simply  rely 
on  computational  verification  that  the  nodal  sets  indeed  allow  for  the  computation  of  a  unique  interpolation 
polynomial  [29].  Under  this  assumption  we  can  likewise  express  Eq.(12)  as 

N 

(13)  Vi  :  ./'(C)  , 

3=0 

which  has  to  be  true  for  any  /  6  C [I] ,  and  in  particular  p,(C  itself.  Hence,  the  Lagrange  polynomials  can  be 
evaluated  at  any  point,  £  €  I,  by  solving  the  dual  problem 


(14) 


\T  L  =  p 
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where  L  =  [Lo(£), Ln(£)]t  and  p  =  [po(£),  ■■,Pn(£)]T ■  This  naturally  enables  the  evaluation  of 
anywhere  in  I  by  computing  Lj{£)  and  applying  Eq.(13). 

In  seeking  the  approximate  solution  to  partial  differential  equations,  the  single  most  important  operation 
is  that  of  computing  approximations  to  spatial  derivatives.  However,  once  we  have  identified  a  well  behaved 
Lagrange  basis,  approximations  to  spatial  derivatives  evaluated  at  the  grid  points,  £,:,  is  obtained  directly 
through  matrix-vector  products  as 


d£ 

dZ 


D  «/ 


3  5/ 


n~6 _ £_ 

”  dp 


D 


3d  f 


'TO  _ £_ 

nd( 


=  DC/  , 


where  the  entries  of  the  quadratic  differentiation  matrices  are  obtained  as 


=  dLjM  n"  =  nC  =  dL Mi) 

«  0<  ’  ij  Oil  ’  d( 


The  entries  can  be  computed  directly  by  using  Eq.(14)  and  the  uniqueness  of  the  polynomials  as 


D«  =  P«V_1  ,  D17  =  P^V-1  ,  Dc  =  P^'V1  , 
where  the  entries  of  P^’  vX)  are 


/i^  p€  =  po  =  dPjiSi)  pC  = 

1  (li'  *  Oil  ’  5C 

4.  A  Convergent  Scheme  for  Maxwell’s  Equations.  Having  realized  high-order  formulations  of 
basic  operations  on  the  nodal  tetrahedron,  we  are  now  in  a  position  to  develop  a  scheme  suitable  for  solving 
linear  systems  of  hyperbolic  problems  in  complex  geometries,  exemplified  by  a  scheme  for  solving  Maxwell’s 
equations. 

To  simplify  matters,  let  us  express  Maxwell’s  equations,  Eq.(3),  in  conservation  form 


(16)  -£  +  \7-F(q)  =  S  , 

where  we  have  introduced  the  state  vector,  q,  and  F(q)  =  [Fi(q) ,  F2(q) ,  F3(q)]T ,  as 


the  flux  defined  as 


srE 

prH 


Fdq) 


— e,  x  H 
e,:  x  E 


respectively.  Here  e,  signifies  the  three  Cartesian  unit  vectors  and  S  =  [SE,SH]T  represents  body  forces, 
e.g.,  currents,  and  terms  introduced  by  the  scattered  field  formulation,  Eqs.  (7)-(8). 


4.1.  Central  Elements  of  the  Scheme.  Let  us  begin  by  introducing  the  nodal  basis  discussed  in  the 
previous  section  and  assume  that  the  statevector,  q,  can  be  represented  as 

N 

qN(x,t)  =  '^2q(xj,t)Lj(x)  , 
i=o 

within  each  general  curvilinear  element,  DC 

We  shall  consider  schemes  in  which  we  require  Eq.(16)  to  be  satisfied  in  the  following  way 


(17) 


dq 


N 


dt 


+  V  •  Fn  —  Sn  cpi(x)  dx 


yji(x)C([qN])dx  . 
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Here  fa  and  ipi  signify  sequences  of  N  functions  while  G([qAr])  is  a  function  of  the  jump  [qw]  of  the  statevector 
at  the  boundary /interface  of  the  element,  e.g.,  if  the  face  is  at  a  solid  boundary  [qAr]  reflects  the  difference 
between  the  prescribed  boundary  condition  and  the  actual  value  of  the  statevector. 

Let  us  emphasizing  a  few  characteristics  of  this  general  formulation,  Eq.(17).  In  particular  we  see  that 
consistency  of  the  scheme  is  immediate  as  the  right  hand  side  of  Eq.(17)  vanishes  when  the  exact  solution  is 
introduced,  i.e.,  if  the  inner  scheme  is  consistent  so  is  the  full  approximation.  One  should  also  observe  that 
boundary /interface  conditions  are  not  imposed  exactly  but  rather  weakly  through  the  penalizing  surface 
integral.  Finally  we  emphasize  that  in  a  multi-element  context,  the  formulation  is  inherently  discontinuous, 
enforcing  the  interface  conditions  weakly  through  the  penalizing  term  and  giving  rise  to  a  highly  parallel 
formulation  of  the  scheme. 

In  choosing  <f>i,  fa  and  Gjfqjy])  one  has  a  tremendous  degree  of  freedom  in  designing  schemes  suitable  for 
solving  differential  equations.  In  [10]  we  proposed  stable  spectral  collocation  methods  with  weakly  imposed 
boundary /interface  conditions  for  solving  the  advection-diffusion  equation  and  the  compressible  Navier- 
Stokes  equations  by  choosing  fa(x)  =  tpi(x)  =  S(x  —  Xj)  and  defining  G([qAr])  to  impose  the  correct  upwind 
flux  conditions.  Alternative  choices,  likewise  leading  to  stable  schemes  for  solving  linear  conservation  laws, 
were  discussed  in  [28,  29].  There  we  considered  mixed  Galerkin-collocation  formulations  by  choosing  fa(x)  = 
Lj(x ),  as  in  a  classic  Galerkin  formulation,  but  using  ipi{x)  =  S(x  —  Xj)  to  impose  the  boundary/interface 
conditions.  Upwind  flux  conditions  were  used  to  construct  GQqjy]). 

To  formulate  a  scheme  for  Maxwell’s  equations,  let  us  assume  that  the  electric,  E,  and  magnetic,  H , 
field  components  can  be  represented  as 

N  N 

Ejy(x,ty  =  ^  ]  E(Xj,t)Lj(x)  =  ^  ]  Ej(t)Lj(x)  , 


j= 0 

j= o 

N 

N 

H N(x,t)  =  ^ ~^H(xht)Lj{x )  =  H i(t)Lj(x)  , 

./■•<)  j— o 

within  each  general  curvilinear  element,  Dfc.  Here  Ej(t )  and  Hj(t)  represent  the  time  dependent  nodal 
values,  i.e.,  the  unknowns  of  the  scheme,  while  Xj  =  Xj{£j)  are  the  mapped  nodal  coordinates. 

We  shall  require  that  the  equations,  Eq.(3),  be  satisfied  in  the  following  Galerkin-like  way 

(18)  +  V-Fn-Sn^J  L,{x)  dx  =  f  T(x)Li(x)n  ■  [F+]  dx  , 

where  qN ,  Fjy,  and  Sn  refers  to  the  approximate  state  vector,  flux  and  body  force,  respectively.  As  in  Sec. 
3,  Lj(x)  represents  the  n’th  order  Lagrange  interpolation  polynomial,  i.e.,  in  the  language  of  the  general 
formulation  in  Eq.(17)  we  have  <pj(x)  =  ipi(x)  =  Lj(x),  while  we  have  G([qjv])  =  r(x)n  ■  [J’jy].  Here  n  is 
the  outward  pointing  normal  vector,  t(x)  is  a  free  parameter  to  be  specified  later,  while  [ F jy]  reflects  the 
jump  in  the  upwind  flux,  i.e.,  we  have  introduced  the  splitting,  iTy  =  F ^  +  F^,  into  the  upwind,  FAr,  and 
downwind,  F] component  of  the  flux. 

It  is  noteworthy  that  the  classical  discontinuous  Galerkin  formulation  [23]  is  recovered  from  Eq.(18)  by 
a  simple  integration  by  parts  and  considering  all  fluxes  at  interfaces  as  upwind  fluxes,  i.e.,  it  is  a  special  case 
of  the  much  more  general  approach  put  forward  in  Eq.(17). 

To  understand  the  exact  form  of  the  penalizing  flux  term,  n  ■  [f  Ar],  it  is  helpful  to  recall  that 


11 


n-  FN  = 


—n  x  Hn 


n  x  EN  | 

i.e.,  the  normal  component  of  the  fluxes  represents  nothing  else  than  the  tangential  field  components  and 
the  effect  of  the  right  hand  side  in  Eq.(18)  is  to  impose  the  correct  boundary /interface  conditions  on  the 
tangential  field  components  at  the  face  of  the  element.  It  is  worth  noticing  that  the  unspecified  function, 
t(x),  controls  how  strongly  the  conditions  are  enforced,  e.g.  if  r  is  very  large  the  conditions  are  essentially 
enforced  exactly. 

As  discussed  in  Sec.  2  the  boundary  conditions  on  the  tangential  field  components,  be  that  in  the 
scattered  field  or  in  the  total  field  formulation,  require  continuity  between  any  two  elements  regardless  of 
their  material  properties.  This  yields  the  explicit  form  of  the  penalizing  boundary  term  as  [32] 


(19) 


where 


n  '  K] 


Z  1h  x  ( Z+[Hn ]  -fix  [Ebv]) 
F_1n  x  (— n  x  [Hn]  -  I'+I-Eov]) 


[En]  =  E+  -  E~n  ,  [H N]  =  H+  -  Tf-  , 

measures  the  jump  in  the  field  values  across  an  interface,  i.e.,  superscript  refers  to  field  values  from  the 
neighbor  element  while  superscript  refers  to  field  values  local  to  the  element.  To  account  for  the  potential 
differences  in  material  properties  in  the  two  elements,  we  have  introduced  the  local  impedance,  Z ±,  and 
conductance,  F±,  defined  as 


and  the  sums 


Z  =  Z+  +  Z~  ,  Y  =  Y+  +  Y~  , 

of  the  local  impedance  and  conductance,  respectively. 

The  special  case  of  a  perfectly  conducting  wall  is  handled  in  the  above  formulation  be  defining  a  mirror 
state  within  the  metallic  scatterer  as 

n  x  E^r  =  -n  x  E] ^  ,  n  x  H^,  =  h  x  H ^  , 

to  enforce  the  correct  boundary  conditions  and  define  the  material  parameters  by  Z+  =  Z~ . 

Now  returning  to  the  semi-discrete  scheme,  Eq.(18),  we  have  an  elementwise  expression  for  the  electric 

field 


N 


(20) 


£  ( M|  dE‘ 


j= 0 


dt 


~  S ij  x  H j  -  MijSf) 


=  EF« 


a  n,  x 


Z+[H,]  -n,  x  [E, 

z^  +  zr 
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and  likewise  for  the  magnetic  field  components 


(21) 


j—0 


E(M«^r  +  s«x^-M«sf) 

-ni  x  [H,]  -  Y+jE,' 


=  E  Fi>  ( 


E  + 


Here  we  have  introduced 


Mf )  =  (Li(x),e(x)Lj(x))D  ,  M'f  =  {Ll{x)^{x)Lj{x))D  , 

as  the  material  scaled  mass-matrices  and 

M  ij  =  (L;(x),Lj(x))d  ,  S;;  =  ( Si ; .  S", .  S )  =  (L,;  (*) ,  VLj  (*))  D  , 

representing  the  local  mass-  and  stiffness  matrix.  Note  that  in  the  special  case  where  er  and  pT  are  elemen¬ 
twise  constant,  we  recover  (XT'.  VI'1)  =  (erM,pTM). 

We  have,  furthermore,  introduced  the  face-based  mass  matrices 


F  u  =  (Li(x),T(x)Li(x))SD  , 

where  the  second  index  is  limited  to  the  trace  of  the  nodal  set  situated  at  the  faces  of  D. 
Expressing  Eqs.(20)-(21)  in  fully  explicit  form  yields 


(22) 


and 


dE 


N 


=(Me)_1  S  x  Hn  +  (M£)_1  MS 
+  (Mf)_1  F  (nx  Z+iHN]-nx[EN] 


-l 


dt 


Z+  +  z- 


rtD 


(23) 


dH  N 
dt 


(M")_1  S  x  En  +  (M^)-1  M Sh 


(M'‘)_1F 


/  nx[HN]  +  Y+[EN] 

n  x  - - — H - - - - 

V  Y+  +  Y- 


<5D 


The  discrete  operators  that  need  to  be  initialized  are,  besides  the  mass-matrices,  M  and  M^,  which  can  be 
computed  exactly  as  described  in  the  appendix  and  inverted  straightforwardly.  We  shall  also  need 

(M-'T1  S  =  (M^'T1  [S'r,  S»,  S:]T  , 


representing  the  general  curvilinear  differentiation  matrix,  as  well  as  (M^)-1  M  for  source  terms.  It  is 
worth  noticing  that  for  all  straightfaced  tetrahedra  with  constant  material  parameters,  the  entries  of  S  can 
be  formed  directly  by  combinations  of  the  classical  differentiation  matrices  introduced  in  Sec.  3,  e.g., 


M  'Sf  Dicy  +  |)S/j  +  |),c  s 

and  similarly  for  M_1Sy  and  M_1SS.  Hence,  as  discussed  in  detail  in  the  appendix,  template  matrices  can  be 
used  for  the  initialization  of  these  operators  in  all  such  elements  while  an  individual  initialization  is  required 
for  general  curved  elements  and  elements  with  smoothly  varying  material  parameters. 

The  same  holds  true  for  the  face-based  operators  M-1F  which  again  can  be  precomputed  for  all  straight- 
faced  elements  with  constant  materials  by  linear  scaling  from  standard  template  operators  for  I .  The  general 
curvilinear  faces  requires  individual  attention. 
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4.2.  Consistency.  In  analyzing  the  scheme,  Eqs.(22)-(23),  it  is  natural  first  to  consider  the  global 
accuracy,  and  hence  consistency,  and  how  it  depends  on  the  size  of  the  tetrahedra,  i.e.,  its  /i-convergence 
rate,  as  well  as  how  it  scales  with  the  order,  n,  of  the  polynomial  approximation.  To  simplify  matters  we 
shall  assume  throughout  that  all  elements  involved  are  straightfaced,  i.e.,  the  transformation  between  D 
and  I  is  linear.  We  shall  furthermore  assume  that  the  material  parameters,  er  and  \Lr ,  be  constant  on  each 
element,  but  they  can  vary  freely  between  elements.  We  shall  later  briefly  revisit  the  impact  on  the  results 
of  the  analysis  of  relaxing  these  assumptions. 

Let  us  introduce  the  exact  solution,  q  =  [E,H ],  to  Maxwell’s  equations,  Eq.(3),  as  well  as  its  projection, 
VnQ  =  [ VnE ,  VnH]t ,  on  the  space  spanned  by  n-order  polynomials,  i.e.,  VnQ  £  Pfr  Except  in  very  special 
cases  VnQ  will  generally  be  different  from  the  numerical  solution,  qN  =  [En,Hn]T,  which  is  the  exact 
solution  to  the  discrete  problem,  Eqs.(22)-(23). 

Before  we  continue  we  wish  to  note  that  a  subtle  consequence  of  using  a  purely  nodal  basis,  as  opposed  to 
a  modal  basis,  is  the  introduction  of  a  discrete  aliasing  error  in  the  interpolation  of  the  initial  conditions.  One 
could  avoid  this  by  reading  the  nodal  values  of  the  Galerkin  projection  of  the  initial  conditions,  computed  by 
using  a  quadrature  of  sufficiently  high  order.  However,  if  the  initial  conditions  are  smooth  and  well  resolved 
this  discrete  aliasing  error  is  small  and  we  shall  not  discuss  it  further  in  what  follows. 

As  the  global  error  is  bounded  by  the  sum  of  the  local,  element-wise  errors,  it  suffices  to  consider  the 
latter.  Introducing  the  exact  solution,  q  =  [ E,H ],  to  Maxwell’s  equations,  Eq.(3),  into  the  semi-discrete 
approximation,  Eqs.(20)-(21),  immediately  yields 


(e,:,Ts)d  =  (Li,  V  x  H  -  VNV  x  H)d  +  [Li,  SE  -  VnSe^d  , 

(L,:,Tff)D  =  -  (Li,  V  x  E  -  Vn\7  x  E)d  +  [Li,SH  -  VNS ff)D  , 

where  Tq  =  jr^bT^J  signifies  the  truncation  error  associated  with  the  scheme.  Note  in  particular  that 
the  surface  terms  of  Eqs.(20)-(21)  vanish  identically  as  the  exact  solution  always  has  smooth  tangential 
components  as  dictated  by  the  physics. 

To  bound  the  truncation  error  we  shall  need  the  following  result  [33,  24,  25] 

Lemma  4.1.  Assume  that  u  £  Wp( D),  p  >  0.  Then  there  exists  a  constant ,  C ,  dependent  on  p  and  the 
angle  condition  of  D,  but  independent  of  u,  h  =  diam( D)  ,  and  n,  such  that 


ha~q 

I|u  -  VNu\\Wq(D)  <  C— —  |HIwp(D)  , 

where  a  =  min(p,  n  +  1)  and  0  <  q  <  a. 

Here  we  have  introduced  the  standard  Sobolev  norm 


iMiw?(D) _  y 

|a|<p 


Qot  i  ga  2  gas 


dxai  8xa2  dxaa 


with  the  multi-index,  a  =  (ai, 02,013). 

With  this  result  and  the  use  of  the  Cauchy-Schwarz  inequality  we  immediately  recover  the  consistency 
result 

Theorem  4.2.  Assume  that  the  exact  solution,  q  =  [E,H]T  G  Wp( D),  p  >  1  and  that  the  body  forces, 
Sq  =  .S^J  £  Wp( D);  p  >  0.  Then  there  exists  a  constant,  C,  dependent  on  p  and  the  angle  condition 
of  D,  but  independent  of  q,  h  =  diam(  D),  and  n,  such  that 
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IIt9IId  ^  c  (^nP- 1  ll9llw'p(D)  +  nP  ll^llwpfD)^  ’ 

where  a  =  min (p,n  4- 1). 

Hence,  if  the  solution  is  locally  smooth  we  can  expect  very  rapid  convergence  in  the  order  of  the  approx¬ 
imation  as  well  as  by  decreasing  the  element  size.  In  particular,  if  the  solution  is  analytic  we  can  expect  to 
recover  full  spectral  convergence  provided  the  scheme  is  stable. 

4.3.  Stability.  Let  us  attend  to  the  issue  of  semi-discrete  stability  and  define  the  local  energy 

Ek  =  ^f  {p\H\2  +  e\E\2)  dx  , 

4  Jok 

and  the  associated  global  energy,  E  =  Ek . 

Local  elementwise  semi-discrete  stability  is  stated  as  follows 

Lemma  4.3  (Local  Stability).  Assume  that  a  solution  to  Maxwell’s  equations  exists  on  the  domain  D. 
If  the  faces  of  the  element  reside  away  from  a  perfect  conductor ,  stability  of  the  semi-discrete  approximation 
to  Maxwell’s  equations,  Eqs.(22)-(23),  is  guaranteed  provided 

1 

T  -  3  ' 

In  case  one  of  the  faces  coincides  with  a  perfect  conductor,  stability  of  the  semi-discrete  approximation  is 
guaranteed  if 


t  =  1  . 


Proof.  For  local  stability  away  from  metallic  boundaries,  it  suffices  to  consider  the  question  of  stability  for 
homogeneous  boundary  conditions,  i.e.,  E^  =  H ^  =  0.  Consider  Maxwell’s  equations  on  the  semi-discrete 
form,  Eqs.(20)-(21),  multiply  from  the  left  with  (Ej,Hj)  and  sum  over  all  the  nodes  in  D  to  obtain 


and 


1  d_ 

2  dt 


( En ,  £-Eov)d  —  ( E  n  5  V  x  ffjv)  q+  (^En  ,  SE  j 


tE 


N 


<5D 


n  x 


Z+H jv  —  n  x  E 


N 


z+  +  z- 


dx 


Id 
2  dt 


(HN,pHN)D  =  x  En) d  +  (HN,SH^ 


D 


/  „  (  „  1  +En  +  n  x  Hn  . 

+  <b  tHn  •  n  x  - ...  -  — -  I  dx 

JS  D  V  Y  +  * 

Adding  the  two  contributions  and  applying  the  the  divergence  theorem  yields 


4- Ek=((  (1  -T)n  ■  {Hn  x  En)  dx 
dt  Jsd 

+  ®  (  =En  ■  fi  x  n  x  En  +  =HN  ■  hx  h  x  H N  )  dx 

Jsd  \Z  Y 


+  (en,S°)d  +  (Hn,s»)d  . 
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Using  standard  vector  identities  this  simplifies  as 


To  ensure  semi-discrete  stability  it  suffices  to  require  that 


(24)  (1  -  t)H„REn  +  Le^RtREn  +  =H%RtRHn  >  0  , 

Z  1 

where  we  have  introduced  the  rotation  matrix 


R  =  R(n) 


0  -n~  nv 

nz  0  —nx 
-ny  nx  0 


Expressing  the  quadratic  form,  Eq.(24),  as  qJjAqN  with  A  reflecting  Eq.(24),  one  recovers  the  first  two 
eigenvalues  of  A  as  Ago  (A)  =  0  while  the  remaining  are  given  as 


A3, 4  — 


r(l  +  Z)±  \J r2(l  +  Z)2  +  Z  (— 3r2  -  2r  +  1) 

2l  ” 


and 


As, e  — 


r(  1  +  Y)  ±  yr2(l  +  F)2  +  Y  (-3r2  -  2r  +  1) 

w 


Hence,  a  sufficient  condition  for  stability  clearly  is  that  r  >  0  and  — 3r2  —  2r  +  1  <  0,  i.e. , 

1 


T  -  3 


In  case  a  face  resides  at  a  metallic  conductor  we  employ  the  boundary  conditions 


n  x  En  =  —n  x  E^-  ,  n  x  H N  =  n  x  , 
and  Z+  =  Z~  =  Z,  Y+  =  Y~  =  Y. 

Following  the  exact  same  procedure  as  above,  we  recover  the  constraint 

(1  -  t)HtnRE  +  ^-EtnRt RE n  >  0  . 

ZZj 

Computing  the  eigenvalues  of  the  corresponding  quadratic  form  yields  two  pairs  of  the  form 

A 1  0  A2j3  =  ^  ±  ~  \/t2  +  Z2(r  -  l)2  . 

Clearly,  the  only  way  to  guarantee  positivity  of  the  eigenvalues  and  hence  the  quadratic  form  is  to  choose 
r  =  1.  □ 

The  result  on  local,  elementwise  stability,  only  supplies  a  necessary  but  not  sufficient  condition  for 
stability.  To  understand  the  issue  of  global  stability  we  must  also  consider  the  influence  of  the  coupling 
between  the  individual  elements. 
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Lemma  4.4  (Face  Stability).  Assume  that  a  solution  to  the  Maxwell’s  equations  exists  on  a  domain  con¬ 
sisting  of  two  elements  sharing  one  common  face.  Stability  of  the  semi-discrete  approximation  of  Maxwell’s 
equations ,  Eqs.(22)-(23 '),  on  this  domain  is  guaranteed  provided 


T  =  1  . 


Proof.  Consider  Maxwell’s  equations  on  the  semi-discrete  form,  Eqs.(20)-(21).  Multiply  from  the  left 
with  (Ej,Hj)  and  sum  over  all  the  nodes  in  D  to  obtain 

2~fif  {En’£En)d  =  x  ^n)d  ^  )D 

r  Z+[Hn] -n~  x  [En}\ 

+Lt  n '  r  x  — — J dx  • 

and 


1  d_ 

2  d.t. 


{HN,pHN)D  — 


(H„,VxEJj)d  +  (h-,Sh)d 


Addition  of  the  two  contributions,  application  of  the  divergence  theorem  and  standard  vector  identities  yields 


To  understand  the  stability  of  a  common  edge,  it  suffices  to  consider  the  case  where  SE  =  SH 
the  contribution  from  two  edges,  utilizing  that  n~  =  — n+,  yields 


0.  Adding 


-rE  =  (f  (1  —  r)(n  •  Hn  x  EN  -  h  ■  x  E^) 
at  Jsd 

T  _  _  T  _  _ 

+=[EN]-n  xn  x  [En]  +  =[Hn]  •  n  xn  x[Hn]cIx 
Z  Y 

=-  /  (1  -  r)n~  ■  (H%  x  E%  -  H „  x  E^) 

JSD 

-  =  | n~  x  [_£bv]|2  -  =  | n~  x  [HN]\2  dx  . 

z  y 

A  sufficient  condition  for  this  to  be  negative  is 


(1  -  r)  ( ( H~n)  T  R E~n  -  {H+f  RE+  )  + 

^[En]tRtR[En]t  +  ^[jrA-]TRTR[LfA,]T  >  0  . 
Z  Y 
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Inspection  reveals  that  by  defining  q  =  [E^ ,  ,  H^]T ,  the  condition  may  be  expressed  is  given  as  a 

symmetric  quadratic  form,  i.e.,  it  suffices  to  choose  r  such  that  all  eigenvalues  of  A  are  non-negative.  Leaving 
out  the  lengthy  and  purely  algebraic  manipulations,  we  consider  the  resulting  two  sets  of  eigenvalues  of  A 
given  as 


Ai,2  =  0  ,  A3,4  =  =  ±-L\/4t2  +  Z(t  -  l)2 

and 


As, 6  =  =  ±  \/4r2  +  Y(t  —  i)2  . 

Clearly,  the  choice  of  r  =  1  is  the  only  feasible  solution  that  ensures  stability  of  the  upwind  scheme  used  for 
connecting  the  elements.  0 

With  these  results  in  place,  we  can  now  state 

Theorem  4.5  (Global  Stability).  Assume  that  a  unique  solution  to  Maxwell’s  equations  exists  in  the 
general  domain,  0.  Assume  furthermore  that  the  boundary  of  Q  is  either  periodic  or  terminated  with  a 
perfectly  conducting  boundary. 

Then  the  semi-discrete  approximation  to  Maxwell’s  equations,  Eqs.(22)-(23),  is  globally  stable  in  the 
sense  that 


4- E<C 
dt 


E  + 


+ 


-,H 


provided  only  that 


t  =  1  . 


Proof.  As  each  face  is  counted  only  once,  the  result  follows  directly  by  summation  over  the  all  the  faces 
and  the  application  of  Lemma  4.3  and  Lemma  4.4 


dt 


E<Y,(EN.SE)Di  +  (HN.SH)t 


<C  [E  + 


+ 


-,H 


+ 


),  ||£1jv||d  <  C  (EpfisrE]y)D  since  £  >  1.  A  similar  line  of 


using  that  ^.Eov,iS£'j  <  C'dl^jv | 
reasoning  is  applicable  for  ^iTjv,  SH^j  and  the  result  on  global  stability  follows.  □ 

4.4.  Convergence.  Having  established  consistency  as  well  as  stability  in  equivalent  norms,  convergence 
follows  directly  from  the  equivalence  theorem  with  a  bound  on  the  local  error 


eD(f)  —  \\E(t)  -  E N 


+  II  H(t)-H 


mciiD 


of  the  form 


eD(t)  <  Ceat  (eD(0)  +  J*  ||T«(S)||D  ds 
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and  global  convergence  is  hence  established  up  to  exponential  growth  in  time  as  is  typical  for  Lax-type 
stability  results. 

As  it  turns  out,  however,  we  can  do  better  and  recover  a  sharp  bound  for  the  growth  in  time  by 
generalizing  ideas  first  put  forward  in  the  context  of  finite  difference  methods  [34],  To  realize  this,  let  us 
make  the  natural  split  of  the  elementwise  error  as 


eD  <  (|| E-  VnE\\d  +  |1  H  -  VnH ||d)  +  (\\VNE 


=  cn  +  £b 


D  > 


EN\\D  +  IjiPjv-ZV  -  H a7 1 1 D ) 


where  cq  is  due  to  the  error  introduced  by  the  polynomial  approximation  of  the  exact  solution  while  cq 
measures  the  errors  associated  with  the  semi-discrete  approximation  of  Maxwell’s  equations. 

To  bound  we  need  only  recall  Lemma  4.1  to  state 

Lemma  4.6.  Assume  that  q  =  [E,H]T  G  WP( D).  Then  there  exists  a  constant,  C ,  dependent  on  p  and 
the  angle  condition  of  D,  but  independent  of  q,  h  =  diam( D),  and  n,  such  that 

ha 

\\q  -  VnqWd  <  C— \\q\\WP^  > 

where  a  =  min(p,  n  -f  1)  and  p  >  0. 

To  arrive  at  a  bound  for  e^,  let  us  first  consider  the  projection  of  the  truncation  error,  VnTq  = 
\VNTE ,VnTh1\  ,  on  the  form 


(25)  (Li,VNTE^  d  =(Lit  Vn\7  xH-  Vn\7  x  VnH)d 

{Li,n  x  (Z+[PnH\  -  fix  [VnE]))sd  , 

Zj 

(26)  ( LhVNTH^D  =- (Li,VNX7  xE-  VNV  x  VNE)D 

~=  {Li, fix  (-Y+[PnE]  -fix  [PnH]))sd  . 

This  is  derived  by  introducing  VnQ  into  the  semi-discrete  scheme,  Eqs.(20)-(21),  exploiting  that  q  satisfies 
Maxwell’s  equations,  Eq.(3). 

The  projection  of  the  truncation  error  can  be  bounded  by  the  exact  solution  as 

Lemma  4.7.  Assume  that  q  =  [E,H]T  G  Wp(D),p  >  3/2.  Then  there  exists  a  constant,  C,  dependent 
on  p,  the  angle  condition  of  D  and  the  local  material  properties,  er,pr%  but  independent  of  q,  h  =  diam( D), 
and  n,  such  that 

H'PjVT1D  <  C  np-3/2  1 1  ^  1 1 W?  ( D )  > 

where  a  =  min(p,  n  +  1). 

Proof.  We  need  only  establish  the  result  for  VnTe ,  Eq.(25),  as  the  derivation  of  the  result  for  VnTe 
following  identical  lines. 

As  VnTe  G  Pf,  =  J2jTE Lj(x)  we  can  multiply  from  the  left  with  TE  and  sum  over  all  the  nodes  to 
recover 
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VnT 


E 


--(vnTe,VnV  x  (H-VnH 
-=  (vNTE,nx  ( Z+[PnH }  -  n  x  [TNE])) 


5D 


Using  the  Cauchy-Schwarz  inequality  and  the  estimate  [25] 

IIQjvIUd  <  ll^jvll 


D  ' 


for  all  q  G  Pf,  (D),  h  =  diain(D),  we  recover 


(27) 


VnTe 


<Ct\\VNVx  (H-VnH)\\d 

+C^=WZ*\'P»Hr]-[V«Er\\\sa 


where  we  for  simplicity  have  introduced  the  tangential  components 


Et  =  n  x  E  ,  Ht  =  n  x  H  . 
To  bound  the  first  term  we  invoke  Lemma  4.1  to  obtain 


(28)  \\VNV  x(H-  VnH) |]d  <  ||V  x  (H  -  VNH) ||D  <  C ^  \\H\\Wp(D)  . 

Consider  now  terms  of  the  type 


\\[VNET}\\dD  <  \\VnE+  -  E 


T  IUd 


+  \\VnE  -  E 


t  IUd 


where  E +  =  ET  =  ET  represents  the  exact  solution  at  <5D.  Recalling  the  trace  inequality  [35] 

IklliD  <  c  (|MId  llVgllo  +  1  Ikllb)  i  1  e  kC1(D)  , 


implies  that 


\\q -VNqWsD  <  C  (\\q -VNqWoWq -VNqWwno')  +  h  1  II?  -  VNq\\b  )  , 


and  we  recover  by  combination  with  Lemma  4.1  the  bound 

\\[VNET]\\SD  <  Cnp_1/2  \\E\\Wp(D) 
Combining  this  with  Eqs.(27)-(28)  one  obtains  the  result 


VnT 


,E 


\E\\wp(D)  "b  1 1 H 1 1  ppy  { o ) 


)  • 


where  (Ci,Co)  are  independent  of  h  and  n  but  Cb  depends  on  the  local  material  properties  (Z±,Y±). 
The  result  for  ^  rriH^  - - J  n -  - ,J- 


VnT 


D 


is  recovered  in  the  same  way,  yielding  the  result 


VnT 


,  H 


7  —  1 


ha~ 1  ha~ 

D  -  Cl  nP-1  ^E^w*(D)  +  C2  —  3/2 


E\\wr(D)  +  1 1 H 1 1  Wp  { D ) 


)  • 
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hence  establishing  the  stated  result.  0 

Let  us  now  return  to  the  original  quest  for  an  improved  convergence  estimate  and  consider  the  error 
equation 


(29)  ^Li,e—(VNE  -  EN)^j  =(L$^PnX7  x  (VnE  -  En))d 

+  =  ( Li,n  x  ( Z+[PnH  -  Hn]  -  ft  x  [' PNE  -  EN])) SD 

Zj 

p(luvnte)d  , 

for  the  electric  field  and  similarly  for  the  magnetic  field 

(30)  (^4  (Pjv-ff  -  HN)j  =-  (L„VnV  x  (VnH  -  HN))D 

—  =  {Li,  n  x  ( Y+[PnE  -  En]  +  ft  x  [PnH  -  HN]))  SD 
+  (L.l,VNTH^jD  . 

The  combination  of  these  expressions  with  Lemma  4.7  and  the  methodology  of  the  stability  proof  in  Sec. 
4.3  yields  the  improved  convergence  result 

Theorem  4.8.  Assume  that  a  solution,  q  £  Wp( D),  p  >  3/2  to  Maxwell’s  equations  in  fl  =  |Jfc 
exists.  Then  the  numerical  solution ,  qN,  to  the  semi-discrete  approximation  Eqs.(22)-(23)  converges  to  the 
exact  solution  and  the  global  error,  ||q  —  Qjv Hd*  bounded  as 


-qN(t)\\Dk  {WM)  -  vNq(t) ||Dfe 

k  k 


+  \\VNq(0)  ~  Qjv(0)IID^  +*  max  |]T«(s) 

*€[0,i] 


Id* 


<^E 


nP 


|9(0)|| 


WP(Dh 


(Dfe 


where  C  depends  on  the  material  properties  and  the  angle  conditions  of  the  elements  but  not  on  h  and  n. 

Proof.  Since  Vi \E  —  En  &  Pf,  and  VnH  —  Hn  £  Pf,  we  can  use  these  as  elementwise  test  functions  in 
Eq.(29)  and  Eq.(30),  respectively,  to  obtain 


2  ~dt  ^ N ’  ~  En)) d  +  ( VnH  -  HN,  e(VnH  -  Hn))q) 

=  /  [ft  ■  [VnH  -  Hn)  x  (VnE  -  EN) 

J  SD 

+=(VnE  —  En)  ft  x  {ZE [VnH  —  Hn]  —  ft  x  [PnE  —  En]) 

z 

-  =(VnH  -  Hn)  ■  n  X  (Y+[PnE  -EN]+nx  [PNH  -  HN])^j  dx 

(vne-en,te^d  +  (vnh -hn,th^d  , 

where  we  have  employed  integration  by  parts  once.  Following  the  approach  of  Lemma  4.4  we  sum  over  all 
the  faces  to  obtain 
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-  — ^  ((' PnE  -  En,s{VnE  -  EN))Dk  +  (VnH  -  Hn,s(VnH  -  HN))Dk) 

1  k 

<-  E  [^VnE  -  ^JHd‘  +  II [VnH  -  ffjv]|lD‘] 

k 

+  '£[('PNE-EN,VNTE)Dk  +  (vnH-Hn,VnTh}  J  . 

k 

Note  that  since  e  and  p.  are  uniformly  bounded  away  from  zero  the  material  weighted  energy  norm  is  L2- 
equivalent.  Furthermore,  the  term  associated  with  the  jump  at  the  element  interfaces  is  strictly  negative 
and  we  recover  the  bound  on  the  error 

E  W'PnQ  ~  gjvllpfc  <  C  E  ('PnQ  —  Qni  'PnT'1) pt,  , 

k  k 

which,  by  using  the  Cauchy-Schwarz  inequality  and  integration  in  time  yields  the  result 

E  \\vNq(t)  ~  QjvWIId*  <CJ2  fll^ivq(O)  —  Qjv(°)llDfe  +  t  max  ||PjvT9(s)||d^  . 

T  TV  sG[0’i]  ) 

Now  combining  this  with  Lemma  4.6  and  Lemma  4.7  establishes  the  result  and  proves  convergence  on  weak 
assumptions  of  local,  elementwise  smoothness  of  the  solution.  □ 

We  have  hence  established  the  semi-discrete  result  that  the  error  can  not  grow  faster  than  linearly  in 
time  and  that  we  can  control  the  growth  rate  by  increasing  the  resolution.  As  we  shall  verify  in  Sec.  5  this 
linear  growth  is  a  sharp  result.  However,  the  computations  shall  also  verify  that  we  can  expect  that  the 
growth  rate  approaches  zero  spectrally  fast  when  increasing  the  order  of  the  approximation,  n,  provided  the 
solution  is  sufficiently  smooth. 

Prior  to  that,  a  few  comments  are  in  place.  A  rigorous  generalization  of  the  results  obtained  above 
to  cover  situations  with  general  curvilinear  elements  and/or  spatial  variation  of  the  materials  within  each 
element  is  not  straightforward.  This  is  due  to  the  generation  of  higher  order  polynomials  from  the  products 
of  the  individual  polynomial  expressions  of  the  fields,  the  materials  and  the  geometry.  One  can,  however, 
gain  an  intuitive  understanding  of  how  the  geometry  and  material  variations  may  impact  the  accuracy  by 
assuming  that  the  polynomial  representations  are  not  of  the  fields  only  but  rather  of  the  combined  functions, 
\fj(s/e^E,  In  this  case,  we  are  working  only  with  /i-order  polynomial  expansions  and  one  can  expect 

that  the  overall  picture  from  the  results  derived  above  will  hold  approximately  for  these  new  functions. 
Hence,  where  we  originally  had  an  n’th  order  polynomial  to  represent  the  fields,  (E,H)  we  are  now  left 
with  an  n’th  order  polynomial  to  represent  the  combined  variation.  One  consequence  of  this  is  that  we  loose 
accuracy  when  considering  only  the  fields  as  we  essentially  have  to  share  the  resolution  power  between  the 
fields,  the  geometry  as  well  as  the  material  variation.  In  particular,  if  the  element  is  strongly  distorted, 
i.e. ,  J  varies  significantly,  one  can  expect  loss  of  accuracy  as  compared  to  the  straightsided  approximation. 
Provided,  however,  that  the  geometry  is  smooth,  i.e.,  J  nonsingular,  and  the  local  material  variation  is 
smooth,  spectral  convergence  is  preserved. 

4.5.  Convergence  of  Divergence  Error.  In  the  absence  of  sources,  it  is  well  known  that  the  electric 
and  the  magnetic  fields  must  remain  solenoidal  throughout  the  computation.  An  assumption  to  this  effect 
was  indeed  imposed  by  choosing  to  solve  only  Maxwell’s  equations  on  the  form  Eq.(3)  and  considering  the 
divergence  conditions  as  consistency  conditions  on  the  initial  conditions.  However,  given  that  we  can  not 


22 


expect  to  recover  the  projection  of  the  analytic  solution  but  rather  will  compute  a  different,  albeit  convergent, 
solution  we  need  to  consider  the  divergence  of  this  numerical  solution  to  justify  the  original  choice  of  solving 
Eq.(3)  only. 

Using  the  results  of  Sec.  4.4  we  can  state 

Theorem  4.9.  Assume  that  a  solution ,  q  £  WP(D),  p  >  7/2  to  Maxwell’s  equations  in  ff  =  (J*,  D* 
exists.  Then  there  exists  a  constant,  C ,  dependent  on  p  arid  the  angle  condition  of  Dfc  ,  but  independent 
of  q,  h  =  diain( D),  and  n,  such  that  the  divergence  of  the  numerical  solution,  qN,  to  the  semi-discrete 
approximation  Eqs.(22)-(23)  is  bounded  as 


£11X7.^)110, 


h*-1  , 


+  t 


h*-2 

- V77T  max 

nP~ 7/“  «e[o,i] 


\q(s)\\ 


WP(  Dfc 


where  a  =  min (p,n  T  1)  and  p  >  0. 

Proof.  Considering  the  local  divergence  of  H  on  any  D  we  have 


IIV  •  (H  -  HN) ||d  <  || V  •  ( H  -  VNH) ||d  +  || V  •  (' VNH  -  HN) ||D 

The  first  term  we  can  bound  immediately  through  Lemma  4.1  as 

IIV  •  ( H  -  VNH) ||d  <  C ^  \\H\\Wp(D]  , 

where  a  =  min  (p,n  +  1)  and  p  >  1. 

Utilizing  the  inverse  inequality  [25] 
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for  all  Un  G  P2( D),  we  can  bound  the  second  term  as 
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by  combining  the  results  of  Lemma  4.7  and  Theorem  4.8.  An  equivalent  bound  can  be  obtained  for  the 
divergence  of  E jy  in  the  case  of  a  source  free  medium  which,  combined  with  the  above,  yields  the  result.  □ 
As  could  be  expected,  the  result  inherits  the  temporal  linear  growth  from  the  convergence  result  and 
confirms  the  possibility  of  recovering  spectral  convergence  of  the  divergence  under  the  assumption  of  sufficient 
smoothness  of  the  solutions.  It  should  be  noted  that  while  the  result  confirms  high-order  accuracy  and 
convergence,  the  estimate  for  the  actual  convergence  rate  is  almost  certainly  suboptimal  and  leaves  room  for 
improvement. 


4.6.  Entr’acte  on  the  Scattered  Field  Formulation.  Let  us  briefly  return  to  an  analysis  of  the 
scattered  field  formulation  discussed  in  Sec.  2.1,  with  the  modified  scattered  field  equations  given  in  Eqs.(7)- 
(8).  We  recall  that  we  split  the  solution,  q,  as 


q  =  qs+q'  , 
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and  exploit  the  linearity  of  Maxwell’s  equations  to  solve  for  the  scattered  field,  q 9 ,  subject  to  the  forcing  by 
the  incident  field,  ql .  As  discussed  in  Sec.  2.1,  this  does  not  alter  the  scheme  in  any  significant  way  except 
at  metallic  boundaries  where  the  boundary  condition  on  the  electric  field  component  takes  the  form 

n  x  Es^+  =  -n  x  E’fT  -  2'PnE1  , 

in  the  notation  of  Lemma  4.4,  while  the  boundary  condition  on  the  magnetic  field  remains 

n  x  Hs^+  =  n  x  H sfil~  . 

Since  this  constitutes  the  only  difference,  we  can  restrict  the  subsequent  analysis  to  the  case  of  a  metallic 
object  in  vacuum  without  loss  of  generality  as  all  other  complications  are  covered  by  the  analysis  of  the  total 
field  scheme. 

It  suffices  to  consider  the  behavior  of  the  computed  solution  which  can  be  bounded  as  stated  in  the 
following. 

Theorem  4.10.  Assume  that  a  scattered  field  solution,  q8  G  Wp( D)„  p  >  3/2  to  Maxwell’s  equations 
in  fl  =  [Jfc  Dfc  exists,  and  that  the  incident  field  ql  G  Wp( D),  p  >  3/2.  Then  the  energy  of  the  numerical 
scattered  field  solution,  qsN,  to  the  semi-discrete  approximation  of  Eqs.(7)-(7)  is  bounded  as 
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where  C  depends  on  the  material  properties  and  the  angle  conditions  of  the  elements  but  not  on  h  and  n. 

Proof.  The  proof  proceeds  in  a  way  very  similar  to  that  of  Theorem  4.8.  Combining  the  equation  for  the 
scattered  field  solution,  qsN,  with  the  equation  describing  the  projection  of  the  incident  field,  V^q1 \  summing 
over  all  the  faces  and  using  qsN  +  VnQ*  as  the  test  function  we  recover 
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where  the  dissipative  terms  are  gathered  over  the  interior  and  PEC  faces  separately  due  to  different  boundary 
conditions,  while  the  global  sum  involves  the  truncation  error,  VnT9’7  ,  associated  with  the  projection  of  the 
incident  field. 

This  latter  term  can  be  bounded  as  in  Lemma  4.7 
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where  a  =  min  (p,n  +  1). 

Proceeding  as  for  Theorem  4.8  we  subsequently  recover 
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from  which 
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thus  establishing  the  result.  □ 

Hence,  also  the  scattered  fields  remain  bounded  up  to  linear  growth  in  time.  An  interesting  difference 
between  this  result  on  that  of  Theorem  4.8  for  the  total  field  formulation  is  that  the  accuracy  and  growth  rate 
of  the  former  is  controlled  solely  by  the  smoothness  of  the  incident  field  with  the  potential  for  exponential 
convergence  for  sufficiently  smooth  illuminating  fields. 

5.  Validation  and  Performance  of  the  Scheme.  Having  developed  the  complete  formulation  for 
the  time-domain  solution  of  Maxwell  equations,  supported  by  a  thorough  convergence  analysis,  it  is  now 
time  to  consider  the  actual  performance  of  the  computational  framework. 

In  the  following  we  shall  discuss  the  validity  of  the  main  theoretical  results  through  a  few  examples 
as  well  as  exemplify  the  versatility  and  overall  accuracy  and  performance  of  the  complete  framework  for  a 
number  of  benchmarks.  Temporal  integration  of  the  semi-discrete  approximation  given  in  Eqs.(20)-(21)  is 
done  using  a  4th  order,  5  stage  low-storage  Runge-Kutta  scheme  [36]  and  a  stability  limited  time-step  scaling 
as 


At  <  CFLinin  \/£rHr\x\  1  > 

with  y/£rnr  reflecting  the  modified  local  speed  of  light  due  to  materials  and 

|VC|  |Vr?|  |VC| 

*  +  An  +  AC  • 

Here  |  •  |  refers  to  the  absolute  value  of  each  and  of  the  vector  components,  i.e.,  |V£|  =  [|Cail£yU£:|]T- 
Hence,  x  provides  a  measure  of  the  local  grid-distortion  as  a  consequence  of  the  mapping,  Tr,  of  I  into  D,  and 
(AC,  A AC)  measures  axial  distance  separating  neighboring  nodal  points  in  I.  In  this  setting  CFL  typically 
takes  values  of  0(1)  while  the  time  step,  At,  scales  as  At  ~  l/n 2  where  l  is  the  minimum  edge  length  on  all 
tetrahedra  and  n  is  the  polynomial  order  of  the  approximation. 

As  a  general  measure  of  error  we  shall  use  the  discrete  Lp-norm  of  the  error  defined  as 

¥f(t)  1|=  ^E  [fN&j,*) 

where  /jv(a;,t)  is  the  numerical  approximation  to  the  exact  value,  }(x,t)  summed  over  all  nodes,  j,  within 
each  of  the  k  elements. 
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b) 


Fig.  3.  In  a)  is  shown  the  temporal  envelope  of  the  maximum  error  on  Hy(t)  in  the  two-dimensional  cavity  for  different 
orders,  n,  of  the  approximation.  The  slope  of  the  linear  growth  is  plotted  in  b ),  confirming  spectral  convergence  as  predicted 
in  Theorem  4-8. 

5.1.  Elementary  Tests  and  Verification  of  Theoretical  Results.  As  a  first  verification  of  the 
theoretical  estimates,  and  in  particular  the  linear  growth  predicted  in  Theorem  4.8,  we  consider  the  solution 
of  the  two-dimensional  Maxwell’s  equations  in  the  TM-polarization,  i.e. ,  we  solve  for  (Hx,Hy,Ez).  There 
is,  however,  nothing  special  about  this  polarization. 

The  computational  problem  is  that  of  a  simple  two-dimensional  vacuum  filled  cavity,  assumed  to  be 
defined  by  (x. y)  £  [—1, 1]  x  [—0.25,0.25],  with  the  walls  at  x  =  ±1  taken  to  be  perfect  electrical  conducts 
while  the  cavity  is  assumed  to  be  periodic  in  the  (/-direction.  The  initial  condition  is  a  simple  oscillatory 
cavity  solution  as 


Hx(x,  y,  0)  =  0  ,  Hy(x,y,0)  =  cos(irx)  ,  Ez(x,  y,  0)  =  0  , 

and  the  computational  domain  is  discretized  by  8  equivalent  isosceles,  each  with  0.5  wavelength  long  sides. 

In  Fig.  3  we  show  the  temporal  envelope  of  the  maximum  error  of  Hy(t ),  computed  using  the  same 
eight  elements  while  increasing  the  order  of  the  approximation.  Following  the  main  result,  Theorem  4.8,  we 
expect  that  the  error  can  grow  at  most  linearly  in  time  and  that  the  growth  rate  should  vanish  spectrally  for 
smooth  solution.  The  results  in  Fig.  3  not  only  confirm  the  validity  of  both  statements  but  also  illustrates 
that  Theorem  4.8  is  sharp,  i.e.,  we  can  not  in  general  guarantee  slower  than  linear  error  growth,  although 
we  can  control  the  growth  rate  by  the  order  of  the  approximation. 

To  further  evaluate  the  performance  of  the  scheme,  let  us  briefly  consider  the  behavior  of  the  divergence 
and  the  ability  of  the  scheme  to  propagate  waves  over  long  distances.  For  this  purpose  we  shall  continue  to 
consider  the  propagation  of  plane  waves  in  simple  rectangular  domains,  tiled  using  isosceles,  each  with  an 
edge  length  of  0.5  wavelength.  In  Fig.  4  we  show  the  global  L2-error  of  the  divergence  of  H  for  a  plane  wave 
propagating  in  a  fully  periodic  domain  being  2  wavelengths  long  and  0.5  wavelength  wide,  tiled  using  only  8 
triangles.  Consistent  with  the  theoretical  result  in  Theorem  4.9  the  scheme  preserves  the  divergence  error  to 
the  order  of  the  scheme,  i.e.,  the  error  vanishes  spectrally  as  we  refine  the  order,  n,  of  the  approximation.  The 
very  notable  even-odd  behavior  in  the  convergence  is  a  consequence  of  the  alignment  with  the  triangulation. 

The  ability  to  propagate  waves  over  very  long  distances  is  likewise  illustrated  in  Fig.  4  where  we  also  show 
the  T2-error  of  the  Hy  component.  Contrary  to  the  small  problems  considered  first,  we  are  here  considering 
a  200  wavelength  long  domain  and  with  the  exact  solution  being  use  to  truncate  the  computational  domain. 
The  domain  is  tiled  using  isosceles  with  an  edge  length  of  0.5  wavelength  and  a  total  of  800  elements.  We 
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Fig.  4.  In  a)  is  shown  the  global  L2 -error  of  the  divergence  of  H  for  a  plane  wave  propagating  in  a  fully  periodic  domain 
as  a  function  of  time  and  order  of  approximation,  n,  confirming  that  the  scheme  conserves  divergence  to  the  order  of  the 
approximation,  i.e.,  it  decays  spectrally  with  increasing  polynomial  order.  The  L2 -error  of  Hy  as  a  function  of  time  and  order 
of  approximation,  n,  in  a  200  wavelength  long  domain  is  shown  in  b),  confirming  the  ability  to  propagate  waves  over  very  long 
periods  of  time  using  only  few  points  per  wavelength. 


Fig.  5.  In  a)  we  illustrate  the  prism  tiled  using  three  high-order  tetrahedra  while  b)  illustrates  the  maximum  of  Hx  for  a 
(y,  z)-polarized  plane  wave  propagation  as  a  function  of  time  and  order  of  the  approximation,  n,  confirming  spectral  convergence 
for  the  three-dimensional  case. 

observe  in  Fig.  4  an  expected  slow  error  growth  until  t  =  200  after  which  it  settles  at  a  maximum  error 
level.  This  level,  however,  decays  spectrally  as  we  increase  the  order,  n,  of  the  approximation.  Using  as  a 
guideline  that  two  edges  span  a  wavelength,  we  see  that  with  7  points  per  wavelength  (two  n  =  3  triangles) 
yields  about  10%  error,  only  9  points  per  wavelength  (two  n  =  4  triangles)  results  in  about  1%  error  while  11 
points  per  wavelength  (two  n  =  5  triangles)  ensures  about  0.1%  error  after  400  periods.  This  is  a  testament 
to  the  advantage  of  using  a  high-order  framework  for  wave  propagation  problems. 

Let  us  finally  consider  a  simple  three-dimensional  test  case  in  which  we  have  tiled  a  straightfaced  prism 
using  three  straightfaced  tetrahedra  as  illustrated  in  Fig.  5.  The  test  is  that  of  a  plane  wave  propagating 
through  the  prism  with  the  exact  solution  being  used  as  the  boundary  conditions.  As  shown  in  Fig.  5  we 
recover  a  rapid  exponential  convergence  as  the  order,  n,  of  the  approximation  is  increased. 

5.2.  Two-Dimensional  Examples.  Having  verified  the  performance  of  the  basic  computational  setup 
as  well  as  the  theoretical  estimates,  let  us  now  consider  problems  of  a  less  simple  and  more  realistic  character. 
This  shall  not  only  allow  us  to  illustrate  more  general  features  of  the  proposed  framework  but  shall  also  be 
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Fig.  6.  In  a)  is  shown  the  finite  element  grid,  consisting  of  85  f  triangles,  used  for  computing  scattering  by  a  perfect 
electrically  conducting  cylinder  of  size  ka  =  15-7T.  A  section  of  the  grid  in  b)  illustrates  the  body  conforming  nature  of  the  grid 
and  the  nodal  grid  supporting  the  high-order  approximation. 


used  to  verify  that  all  the  properties  of  the  high-order  unstructured  grid  approach,  seen  so  convincingly  in 
the  last  section  for  simple  examples,  carry  over  to  the  solution  of  more  realistic  problems. 

We  shall  focus  the  attention  on  problems  described  by  the  two-dimensional  TM-polarized  Maxwell’s 
equations  on  the  form 


(31) 
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subject  to  boundary  conditions  between  two  regions  with  material  parameters,  er  and  fir  ,  for  k  =  1,  2,  as 


n  x  H(1)  =  nx  H(2)  , 

=  E{2 1  . 

Here  =  (H.k\  Hyk\  0)T  and  n  =  (I) 7  represents  a  unit  vector  normal  to  the  interface.  For  the 

case  of  a  perfectly  conducting  metallic  boundary  the  condition  becomes  particularly  simple  as 

Ez=  0  . 

The  computational  domain  is  truncated  with  a  Cartesian  PML  [37]  using  a  quadratic  absorption  profile. 

It  is  worthwhile  emphasizing  that  results  of  equal  quality  and  overall  accuracy  as  the  ones  shown  in  the 
following  for  the  TM-polarized  case  has  been  obtained  for  the  TE-polarized  case. 

As  a  first  example  we  consider  that  of  plane  wave  scattering  by  a  perfectly  conducting  circular  cylinder 
with  a  radius  of  a  =  7.5A,  i.e. ,  ka  =  157T.  The  surrounding  medium  is  assumed  to  be  vacuum,  i.e. ,  er  =  \ir  =  1. 
The  finite  element  grid,  consisting  of  854  triangles,  utilized  for  this  computation  is  shown  in  Fig.  6  along  with 
a  section  of  the  grid  illustrating  the  full  bodyconforming  nature  of  the  approximation  as  well  as  the  nodal  grid 
supporting  the  high-order  approximation.  Maxwell’s  equations  are  solved  in  the  scattered  field  formulation 
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a)  b) 


Fig.  7.  In  a)  is  shown  the  bistatic  radar  cross  section,  RCS(O),  as  computed  using  the  exact  series  representation  as  well 
as  the  unstructured  grid  method  at  different  polynomial  orders,  n.  Evidence  of  high-order  convergence  for  the  RCS- computation 
is  given  in  b)  showing  exponential  decay  of  the  error  in  RCS  (dBm)  with  increasing  order  of  the  approximation. 


and  Pronv  extrapolation  [38]  is  used  to  reduce  the  required  computing  time  to  reach  the  harmonic  steady 
state. 

In  Fig.  7  we  compare  the  computed  bistatic  radar  cross  section,  RCS(0),  with  the  exact  series  solution 
[39],  for  various  orders,  n,  of  the  approximation  using  the  finite  element  grid  illustrated  in  Fig.  6.  As 
expected  we  observe  a  very  rapid  convergence  with  increasing  n,  yielding  a  reasonable  engineering  accuracy 
computation  with  the  4th  order  scheme  while  increasing  the  order  to  n  =  8  results  in  a  perfect  match.  A 
quantitative  confirmation  of  this  is  also  shown  in  Fig.  7,  illustrating  the  expected  exponential  convergence 
of  the  RCS  with  increasing  n. 

One  of  the  most  appealing  advantages  of  a  high-order  framework  on  simplices  is  the  ability  to  import 
a  strongly  skewed  finite  element  grid  and  recover  a  fully  converged  solution  by  increasing  the  order  of  the 
approximation  rather  than  having  to  reconstruct  an  improved  finite  element  discretization.  This  property 
is  particularly  important  and  useful  for  large  three-dimensional  problems  where  the  grid  generation  phase 
can  be  very  complex  and  time-consuming.  As  an  illustration  of  this  approach  to  convergence,  we  consider 
in  Fig.  8  the  plane  wave  scattering  from  a  PEC  cylinder  with  a  radius  of  one  wavelength,  i.e. ,  ka  =  ‘2ir.  The 
measure  of  accuracy  and  convergence  is  based  on  the  observation  that  the  symmetry  of  the  problem  makes 
one  expect  the  scattered  fields  themselves  maintain  a  high  degree  symmetry. 

This  is  indeed  confirmed  in  Fig.  8  where  we  show  a  deliberately  chosen  poor  grid  and  the  rapid  recovery 
of  the  symmetry  of  one  of  the  scattered  field  components,  Hx,  as  the  order,  n,  of  the  approximation  is 
increased  without  modifying  the  underlying  finite  element  grid.  The  detail  to  which  the  symmetry  is  restored 
is  particularly  noteworthy. 

As  an  illustration  of  the  capability  to  handle  materials  let  us  consider  plane  wave  scattering  by  a  pene¬ 
trable  circular  cylinder  with  a  radius  of  a  =  3.5A  consisting  of  an  ideal  dielectric  with  er  =  2.0,  i.e.,  similar  to 
that  of  glass.  The  problem  is  again  solved  in  a  pure  scattered  field  formulation  and  the  fully  body-conforming 
finite  element  discretization,  consisting  of  a  total  of  1020  triangles,  is  illustrated  in  Fig.  9.  We  note  that 
the  absorbing  PML  layer,  containing  about  2/3  of  the  total  amount  of  triangles  is  unnecessarily  thick  for 
illustration  only  and  can  be  decreased  without  loss  of  accuracy. 

As  is  likewise  illustrated  in  Fig.  9  we  recover  the  full  bistatic  radar  cross  section,  RCS(#),  with  excellent 
correspondence  to  the  exact  solution  [40]  and  quantitative  agreement  over  a  40  db  dynamic  range. 
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Fig.  8.  Example  of  convergence  by  increasing  the  order  of  the  approximation,  n,  on  a  deliberately  chosen  highly  skewed 
finite  element  grid,  illustrated  in  a).  The  convergence  is  illustrated  in  b)-f)  for  increasing  the  order  from  f'th  order  to  12'th 
order,  showing  a  complete  recovery  of  the  expected  symmetry  of  the  scattered  field  component,  II x  ■ 


5.3.  Three-Dimensional  Examples.  As  a  first  verification  of  the  general  three-dimensional  frame¬ 
work,  let  us  consider  plane  wave  scattering  by  a  ka  =  10  perfectly  conducting  sphere,  the  analytic  solution 
of  which  is  given  by  a  Mie-series  [39] . 

We  use  a  fully  bodyconforming  grid  with  a  total  of  3000  elements,  having  an  average  edge  length  at  the 
sphere  of  4A/5.  Contrary  to  the  two-dimensional  case  where  we  used  a  PML  to  truncate  the  computational 
domain  we  choose  in  the  three-dimensional  case  to  embed  the  sphere  in  a  (20A)3  cube  and  employ  stretching 
of  the  elements  as  one  approaches  the  outer  boundary.  The  grid  is  stretched  such  that  the  average  edge  is 
about  2 A  at  the  outer  boundary.  As  in  the  two-dimensional  case,  all  examples  are  done  using  a  4th  order 
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Fig.  9.  Plane  wave  scattering  by  a  ka  =  7tt  dielectric  circular  cylinder  with  a  relative  permittivity  er  =  2.0.  In  a)  we 
show  the  finite  element  discretization  while  b)  shows  a  comparison  between  the  computed  bistatic  radar  cross  section,  RCS(6), 
obtained  with  a  10'th  order  approximation  and  that  recovered  by  evaluating  the  exact  solution. 


a) 


b) 


Fig.  10.  Plane  wave  scattering  by  a  ka  =  10  metallic  sphere  for  a  fixed  grid  and  increasing  order,  n,  of  the  polynomial 
approximation.  In  a)  we  show  the  convergence  of  RCS(0,0)  for  vertical  polarization  (TM),  while  b)  shows  RCS(6,90)  for 
horizontal  polarization  (TE)  of  the  incident  field. 


low-storage  Runge-Kutta  scheme  to  advance  in  time  and  Prony  extrapolation  to  identify  the  solution. 

In  Fig.  10  we  illustrate  the  convergence  of  the  scheme  with  a  fixed  grid  when  increasing  the  order  of 
the  approximation  within  each  tetrahedron.  Even  for  n  =  3,  i.e.,  a  third  order  scheme  with  about  5  points 
per  wavelength,  do  we  compute  a  reasonable  solution  while  increasing  the  order  yields  a  rapidly  converging 
solution  as  one  would  expect. 

As  a  considerably  more  challenging  problem,  let  us  consider  scattering  by  a  perfectly  conducting  business 
card  sized  metallic  plate  as  illustrated  in  Fig.  11.  The  horizontally  polarized  plane  wave  impinges  at  the 
metallic  plate  at  an  almost  grazing  angle,  causing  the  excitation  of  very  strong  waves  along  the  edges  of  the 
metallic  plate.  These  waves  contribute  dramatically  to  the  scattering  process  and  need  to  be  resolved  to 
accurately  predict  the  far  field  scattering. 

This  problem,  being  one  of  the  EMCC  benchmark  problems  [41]  for  code  validation,  is  addressed  by 
using  a  total  of  27000  straightsided  tetrahedra,  each  supporting  a  4th  order  polynomial  approximation.  The 
average  edge  length  at  the  edge  of  the  business  card  is  approximately  A/5.  The  metallic  plate  is  embedded 
in  a  (20A)3  cube,  with  the  elements  being  stretched  to  about  4A  at  the  outer  boundary. 
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Fig.  11.  In  a)  we  show  the  geometry  for  the  plane  wave  scattering  by  a  metallic  business  card  while  b)  shows  the  comparison 
between  monostatic  RCS  experimental  results  [41]  (full  line)  for  horizontal  polarization  of  the  illuminating  field  and  particular 
computed  data  points  (■). 

In  Fig.  11  we  also  show  the  comparison  between  the  experimentally  measured  monostatic  RCS  [41]  and 
a  number  of  particular  computed  data  points.  Again  we  observe  good  agreement  over  the  full  azimuthal 
range  with  results  well  within  the  experimental  error.  The  most  significant  discrepancy  of  a  few  dB  for  <j>  ~  0 
is  consistent  with  other  published  results  [41]. 

As  a  final  example  of  the  performance  of  the  three-dimensional  framework  we  shall  consider  plane  wave 
scattering  from  a  dielectric  cylinder  of  finite  length.  As  illustrated  in  Fig.  12,  the  length  of  the  cylinder  is  5A 
and  the  non-magnetic  material  has  a  permittivity  of  er  =  2.25,  similar  to  that  of  glass.  Clearly,  the  nature 
of  the  fields  is  less  dramatic  than  in  the  previous  case  and  we  find  that  using  a  total  of  approximately  67000 
elements,  supporting  a  4th  order  approximation  and  with  an  average  vacuum  edge  length  at  the  cylinder  of 
A/3,  suffices  to  accurately  predict  the  far  field  scattering.  The  full  computational  domain  is  a  cylinder  of 
radius  16A  and  length  23A  with  the  stretched  elements  having  a  average  length  of  4A  at  the  outer  boundary. 

In  Fig.  12  we  show  a  direct  comparison  between  the  full  bistatic  RCS  for  a  plane  wave  impinging  directly 
at  the  end  of  the  cylinder  as  computed  using  the  current  framework  as  well  as  an  independently  verified 
pseudospectral  multi-domain  axi-symmetric  code  [12],  As  expected  we  find  an  almost  perfect  agreement 
between  the  results  of  the  two  schemes  over  approximately  50  dB  dynamical  range. 

5.4.  Parallel  Performance.  The  discontinuous  element  formulation  of  the  scheme  enables  a  highly 
efficient  implementation  at  contemporary  large  scale  distributed  memory  machines.  While  this  is  a  lesser 
concern  for  the  two-dimensional  schemes,  it  is  essential  to  enable  the  modeling  of  large  scale  three-dimensional 
problems. 

The  developed  schemes  are  implemented  in  a  combination  of  Fortran  and  C  with  all  computationally 
intensive  part  written  in  Fortran  and  taking  advantage  of  Level  3  BLAS  [42]  where  possible.  The  parallel 
interface  is  written  in  MPI  [43]  with  METIS  [44]  used  to  distribute  the  elements  over  the  processors.  To 
ensure  high  cache  efficiency,  we  employ  bandwidth  minimization  [45]  of  the  nodal  points  locally  to  the 
processors  [46].  For  computations  maximizing  the  capacity  of  the  processors,  i.e. ,  filling  the  local  memory, 
this  is  critical  to  ensure  high  performance. 

In  Table  1  we  list  the  parallel  speedup  relative  to  the  n  =  2  case  as  the  number  of  processors  are 
increased.  A  few  things  are  worth  noting.  For  a  fixed  size  problem,  the  parallel  speedup  decreases  slightly 
as  the  number  of  processors  increases  which  is  natural  as  the  relative  communication  cost  increases.  On 
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Table  1 


Parallel  speedup  for  a  123.000  element  grid,  scaled  to  timing  for  n  =  2  on  processors  (-  implies  insufficient  memory 
local  to  the  nodes). 


Polynomial 

Degrees  of 

Number  of  processors 

order  (n) 

freedom  (xlO6) 

4  8 

16 

32 

64 

2 

7.4 

1.0  2.0 

3.9 

7.5 

13.7 

3 

14.8 

0.9 

1.8 

3.5 

6.4 

4 

25.8 

- 

1.0 

1.9 

3.6 

5 

41.3 

- 

- 

0.8 

1.6 

Fig.  12.  In  a)  we  show  the  geometry  for  the  plane  wave  scattering  by  a  dielectric  finite  length  cylinder  while  b)  shows 
the  RCS(O.O)  for  vertical  polarization  (■)  of  the  illuminating  field  and  RCS(0,90)  for  horizontal  polarization  (■)  compared  with 
results  obtained  using  a  pseudospectral  axi-symmetric  code  (full  line)  [12] 

the  other  hand,  for  problem  sizes  utilizing  the  available  resources  we  find  a  very  high  parallel  efficiency, 
e.g.,  increasing  the  problem  size  and  the  number  of  processors  yields  a  close  to  constant  speedup.  The  data 
also  show  a  minor  decrease  in  relative  performance  for  high  order  on  many  processors,  which  we  speculate 
is  related  to  cache  effects  known  to  be  become  important  as  the  size  of  the  operators  increase  [29].  We 
generally  observe  better  than  90%  parallel  efficiency,  consistent  with  other  similar  studies  [47]. 

6.  Concluding  Remarks  and  Outlook.  The  main  purpose  of  paper  has  been  to  introduce  the  reader 
to  a  new  class  of  high  order  unstructured  grid  methods  suitable  for  the  time-domain  solution  of  Maxwell’s 
equations.  A  number  of  central  elements  separate  the  current  framework  from  previous  attempts  to  develop 
high-order  accurate  methods  on  unstructured  grids.  The  use  of  a  purely  nodal  basis  has  a  number  of 
advantages  in  terms  of  ease  of  implementation  by  simple  matrix-vector  operations  as  well  as  the  promise 
to  yield  a  highly  efficient  implementation.  Furthermore,  the  generalized  discontinuous  penalty  scheme  was 
introduced,  offering  an  inherently  parallel  discontinuous  formulation  with  a  purely  block-diagonal  mass 
matrix  which  can  be  inverted  in  preprocessing. 

The  particular  focus  on  Maxwell’s  equations  allowed  us  to  develop  a  complete,  if  not  optimal,  convergence 
theory.  A  similar  analysis  can  be  completed  for  other  classes  of  linear  problems  such  as  acoustics  and  linear 
elasticity.  We  have  confirmed  the  results  of  the  analysis  by  thorough  computational  experiments,  illustrating 
the  flexibility,  versatility,  and  efficiency  of  the  proposed  high-order  accurate  unstructured  grid  framework. 

While  we  have  focused  on  linear  systems  in  general  and  Maxwell’s  equations  in  particular,  the  central 
elements  of  the  framework  allows  for  more  general  formulations  that  enable  the  solution  of  typical  nonlinear 
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systems  of  conservation  laws.  This  naturally  raises  questions  about  proper  formulation  of  the  fluxes  at 
interfaces,  conservation,  entropy  solutions  and  stability  of  high-order  schemes  when  approximating  problems 
with  discontinuous  solutions.  We  shall  address  these  issues  in  [30]  where  we  shall  also  demonstrate  the 
performance  of  such  generalized  formulations  for  the  solution  of  conservation  laws. 

Acknowledgment.  The  authors  extend  their  appreciation  to  Prof.  D.  Gottlieb  and  Dr.  A.  Ditkowski, 
Brown  University,  for  many  fruitful  discussions. 

Efficient  and  Accurate  Implementation  Techniques.  From  the  discussions  in  Sec.  3.2  it  is  clear 
that  the  Vandermonde  matrix,  V,  plays  a  crucial  role  in  setting  up  the  discrete  operators  for  interpolation 
and  differentiation.  The  properties  of  V,  e.g.,  its  conditioning,  depends  exclusively  on  the  structure  of  nodal 
set,  £  •,  and  on  the  way  in  which  we  choose  to  represent  the  basis,  i.e.,  Pi(£).  While  the  former  is  chosen  to 
ensure  well  behaved  Lagrange  interpolation  polynomials,  we  have  significant  freedom  in  the  specification  of 
Piit)- 

A  particularly  simple  choice  is  that  of  the  multivariate  monomial  basis,  i.e.,  Pi(£)  =  .  However, 

even  for  interpolation  in  one  dimension,  i.e.,  Pi(£)  =  £*,  is  it  well  known  that  this  basis  leads  to  the  classical 
Vandermonde  matrix  with  an  exponentially  growing  condition  number.  Hence,  even  for  moderate  values  of 
n  can  we  expect  severe  problems  when  attempting  to  compute  the  action  of  V-1.  The  well  known  solution 
to  this  problem  is  to  choose  a  basis  that  is  orthonormalized  with  respect  to  some  proper  inner  product  to 
assure  the  maximum  degree  of  linear  independence  of  the  basis. 

Such  a  basis  has  been  known  for  long  [48,  49,  50]  and  takes  the  form 


(32) 
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3(2*+2j  +  1,0) 
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where 


20+a_  2(1+1?) 

n+C  ’  i-C 


t  =  C  , 


and  Pna,l3\x)  signifies  the  classical  Jacobi  polynomial  of  order  n  [51]. 

The  tensor  product  structure  of  the  basis,  Eq.(32),  becomes  evident  when  one  realizes  that  while  £  is 
restricted  by  I,  the  mapped  coordinates,  (r, s,t),  covers  [— 1,1]3.  Furthermore,  it  is  easy  to  see  that  the 
polynomial  space  Pf,  can  expressed  as 


pf,,  =  span  {ipiJk  (£) ;  i,  j,  k  >  0j.|  +  j  +  k  <  n}  . 

An  important  property  of  the  basis,  Eq.(32),  is  its  orthogonality  on  I  [21]  as 

Uyfc(£)U  pqr  (4)  ,pqr  5 

where  Sijk,pqr  is  the  multi-dimensional  Dirac  delta  and  the  normalization  is 

2  2  “?+2 

^  J  2i  +  1  2 (i  +  j)  +  2  2 (i  +  j  +  k)  +  3 

Let  us  introduce  the  index,  a  £  [0,iV],  reflecting  some  chosen  ordering  of  ( i,j,k )  and  hence  ipijk-  We  can 
thus  rename  the  polynomial  basis  =  ipa{£)  to  simplify  the  notation  in  the  subsequent  discussion. 
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With  this  machinery  in  place,  let  us  address  how  to  initialize  the  basic  operations  and  the  associated 
operators  needed  for  solving  partial  differential  equations  with  the  current  context  in  an  efficient  and  accurate 
manner. 

Using  the  orthogonal  basis,  ipa,  it  is  natural  to  define  the  Vandermonde  matrix  to  have  the  entries 

V*j  =  ipj(€j)  ■ 
s/Tj 

The  relation  between  the  nodal  and  the  modal  representation  of  a  function,  /,  follows  directly  from  Eq.(12) 
as 


/  =  V/  ,  }  =  . 

Furthermore,  we  can  compute  the  entries  of  the  differentiation  matrices  directly  by  defining  the  entries  of 
p(€,»i,e),  Eq.(15),  using  the  derivatives  of  ?/>,(£)  expressed  explicitly  by  the  identity  [51] 

In  an  equally  simple  and  straightforward  way  we  can  define  spatial  filtering  matrices,  F,  as 

F  =  V<r(i,  j.AOV"1  , 


where  the  order  p  filter  itself  is  defined  as 

a{i,j,k)  =  exp 


(i  +  j  +  k)(i  +  j  +  k  +  3) 


such  that  filtering  is  accomplished  through  a  straightforward  matrix  multiply  at  a  cost  equivalent  to  that  of 
computing  a  spatial  derivative. 

While  the  interpolation,  differentiation,  and  filtering  operators  will  play  a  crucial  role  in  the  solution  of 
the  partial  differential  equations,  we  shall  also  need  to  evaluate  inner  products  on  the  general  curvilinear 
tetrahedron,  i.e.,  we  shall  need  an  efficient  and  accurate  procedure  for  computing 


d  =  J fN(£)gN(€)J{€)d£  , 

where  J  refers  to  the  transformation  Jacobian  for  the  mapping  between  D  and  I  and  /at  G  P%,  9n  G  P„. 
To  evaluate  this  inner  product,  we  exploit  that  /at  and  gjy  are  expressed  uniquely  by  their  expansion  in 
Lagrange  polynomials  as 

N 

(}n,9n) D  =  5Z  frfi 
ij— 0 

Furthermore,  using  the  basis  itself,  wQ  (£),  we  can  express  the  Lagrange  polynomials  themselves  using  Eq.(14) 
on  the  form 


fLd&LjMJMdt 


N 

A;(0  =  EV7^(0 

k= 0 

This  immediately  yields  the  expression 
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(33) 


N  N  r 

( fN,gN)D  =  E  fisj  E  % (tmt)j(t)dt 

i,j= 0  k,l=0 

N  N 

=  E  M'  E  V//Vw"  * 

«,i— 0  k,l=0 

where  the  symmetric  matrix  of  weights,  W,  has  the  entries 

w k,.  =  . 

On  matrix  form  Eq.(33)  becomes 

(/jv,3a0d  =  (V-1/)rWV"1fl  . 

For  all  elements  we  may  precompute  (V_1)TWV_1  in  a  preprocessing  stage,  storing  only  the  upper  half  of 
the  operator  due  to  symmetry.  In  the  particularly  important  case  where  D  is  a  straightsided  tetrahedron, 
i.e. ,  J  is  a  constant,  the  orthonormality  of  ipa  implies  that  W  =  J I,  where  I  represents  the  identity  matrix. 
Hence,  through  a  simple  linear  scaling  one  recovers  the  weights  for  all  tetrahedra  with  planar  faces.  For  the 
general  case  where  J(£)  is  non  constant,  the  entries  of  W  are  computed  exactly  through  over-integration  by¬ 
product  rules  based  on  Legendre  Gauss  quadratures  [52], 

A  final  key  operation  needed  for  the  implementation  of  the  scheme  is  surface  integration,  i.e., 

( fN,9N)sD  =  /  /at(£)<7jv(£)J(0  d£  , 

J  5 1 

where  J(£)  refers  to  the  surface  Jacobian  only.  While  one  could  proceed  as  for  the  volume  integral  discussed 
above,  it  is  more  natural  to  exploit  the  uniqueness  and  completeness  of  the  Lagrange  interpolation.  To 
illustrate  the  procedure,  let  us  restrict  attention  to  one  of  the  faces,  face  ’d’  (see  Fig.  1),  and  term  those 
JVd  =  h(n  +  1  )(n  +  2)  nodes  positioned  at  that  face  for  £d.  Clearly,  using  the  exact  same  procedure  as  for 
the  three-dimensional  Lagrange  polynomial  discussed  above,  we  can  compute  a  two-dimensional  Lagrange 
polynomial,  Zd(£,'/?)  based  on  £d.  As  for  Lj(£),  we  can  recover  Zd  as  the  solution  to  the  dual  problem 

(Vd)T/d  =  pd  , 

where  the  entries  of  the  Vandermonde  matrix  is 


The  proper  basis  to  use  is  the  two-dimensional  version  of  Eq.(32)  given  directly  as  pj(£,  >l)  =  i  '/■  — 1). 
This  allows  us  to  proceed  exactly  as  for  the  volume  integration  and  express  the  integration  over  face ’d’  as 

/  fN(t, n,  -i )gN(t  V,  v,  -1)  dn  =  ( (vd) _1  fd) T  wd  (vd) _1  gd  , 

J  face  d  v  y 

where  fd  =  [/jv(^d),  ■■■,/iv(£ %d)]T  is  tlie  trace  of  Jn  at  the  face.  A  similar  definition  is  used  for  gd.  The 
matrix  of  surface  weights  are  given  as 

wd  =  /  ipiit  V,  -1  )ipjit  V,  -1  )J(Z,  V,  -1)  d^drj  . 

J  face  d 
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In  the  important  special  case  where  the  face  is  planar  and  has  straight  edges,  orthonormality  of  the  poly¬ 
nomials  immediately  implies  that  Wd  =  JdI  as  for  the  volume  case.  For  the  general  case  we  shall  use  a 
cubature  rule  [53,  54,  55]  of  sufficiently  high  order  to  evaluate  the  inner  product,  i.e. ,  we  need  to  interpolate 
the  polynomials,  Jn  and  gw,  onto  the  M  cubature  nodes,  £d,cubj  situated  at  the  face.  This  is  done  by  the 
introduction  of  the  interpolation  operator 

H  I’'  (V1)  '  P  ij=pt(ZfCUb)  , 

i.e.,  P  is  an  Nd  x  M  operator.  The  evaluation  of  the  inner  product  is  then  accomplished  as 
/  fN(Z,ri,-l)gN(Z,ri,-l)J(Z,ri,-l)dSdri=  (fA)TnTWngA  , 

J face  d  V  ' 

where  the  entries  of  the  diagonal  M  x  M  matrix  of  weights  are  given  as 

Nn 

k= 0 

containing  the  weights  w,;  of  the  cubature  as  well  as  the  interpolation  of  the  transformation  Jacobian  of  the 
curvilinear  face.  While  this  formulation  leads  to  the  most  compact  scheme  it  proves  advantageous  to  operate 
directly  on  the  values  at  the  cubature  nodes  as  they  do  not  include  the  edges  and  vertices,  i.e.,  we  can 
establish  a  clean  face  based  connection  between  elements  without  considering  the  multiplicity  of  solutions  at 
vertices  and  the  added  complexity  this  introduces  for  the  implementation  and  performance.  Needless  to  say, 
the  whole  discussion  for  the  evaluation  of  the  integral  over  face  ’d’  carries  over  directly  to  the  other  faces, 
hence  completing  the  evaluation  of  the  full  surface  integral. 

It  is  important  to  realize  that  all  the  operators  introduced  in  the  above  can  be  initialized  during  a 
preprocessing  phase.  Furthermore,  it  is  worth  recalling  the  discussion  in  Sec.  3.1  in  which  we  found  that  any 
two  straightfaced  tetrahedra  are  connected  through  a  linear  transformation.  Hence,  for  any  straightfaced  D 
we  can  form  any  of  the  operators  discussed  in  the  above  directly  by  a  linear  scaling  of  hard-coded  template 
operators  defined  on  I.  This  saves  not  only  preprocessing  time  but  also  reduces  the  required  storage  space 
very  substantially. 
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